]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ANALYSIS/AliESDtrackCuts.cxx
Revert unwanted changes concerning PID in HLT
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDtrackCuts.cxx
index ec144cc15085b4ff79a393d60ba9d068a1845cea..1a07eee202f996a256b3f4ad36c7ffbd478e43c6 100644 (file)
 #include "AliESDtrackCuts.h"
 
 #include <AliESDtrack.h>
-#include <AliESD.h>
+#include <AliESDVertex.h>
 #include <AliESDEvent.h>
 #include <AliLog.h>
 
 #include <TTree.h>
 #include <TCanvas.h>
 #include <TDirectory.h>
+#include <TH2F.h>
+#include <TF1.h>
 
 //____________________________________________________________________
 ClassImp(AliESDtrackCuts)
@@ -35,8 +37,8 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  "require ITS refit",
  "n clusters TPC",
  "n clusters ITS",
- "#Chi^{2}/clusters TPC",
- "#Chi^{2}/clusters ITS",
+ "#Chi^{2}/cluster TPC",
+ "#Chi^{2}/cluster ITS",
  "cov 11",
  "cov 22",
  "cov 33",
@@ -50,8 +52,18 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  "p_{x}",
  "p_{y}",
  "p_{z}",
+ "eta",
  "y",
- "eta"
+ "trk-to-vtx max dca 2D absolute",
+ "trk-to-vtx max dca xy absolute",
+ "trk-to-vtx max dca z absolute",
+ "trk-to-vtx min dca 2D absolute",
+ "trk-to-vtx min dca xy absolute",
+ "trk-to-vtx min dca z absolute",
+ "SPD cluster requirement",
+ "SDD cluster requirement",
+ "SSD cluster requirement",
+ "require ITS stand-alone"
 };
 
 //____________________________________________________________________
@@ -68,8 +80,14 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   fCutAcceptKinkDaughters(0),
   fCutRequireTPCRefit(0),
   fCutRequireITSRefit(0),
+  fCutRequireITSStandAlone(0),
   fCutNsigmaToVertex(0),
   fCutSigmaToVertexRequired(0),
+  fCutMaxDCAToVertexXY(0),
+  fCutMaxDCAToVertexZ(0),
+  fCutMinDCAToVertexXY(0),
+  fCutMinDCAToVertexZ(0),
+  fCutDCAToVertex2D(0),
   fPMin(0),
   fPMax(0),
   fPtMin(0),
@@ -104,9 +122,14 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   SetMaxCovDiagonalElements();                                     
   SetRequireTPCRefit();
   SetRequireITSRefit();
-  SetAcceptKingDaughters();
-  SetMinNsigmaToVertex();
-  SetRequireSigmaToVertex();
+  SetRequireITSStandAlone(kFALSE);
+  SetAcceptKinkDaughters();
+  SetMaxNsigmaToVertex();
+  SetMaxDCAToVertexXY();
+  SetMaxDCAToVertexZ();
+  SetDCAToVertex2D();
+  SetMinDCAToVertexXY();
+  SetMinDCAToVertexZ();
   SetPRange();
   SetPtRange();
   SetPxRange();
@@ -114,6 +137,9 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   SetPzRange();
   SetEtaRange();
   SetRapRange();
+  SetClusterRequirementITS(kSPD);
+  SetClusterRequirementITS(kSDD);
+  SetClusterRequirementITS(kSSD);
 
   SetHistogramsOn();
 }
@@ -132,8 +158,14 @@ AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
   fCutAcceptKinkDaughters(0),
   fCutRequireTPCRefit(0),
   fCutRequireITSRefit(0),
+  fCutRequireITSStandAlone(0),
   fCutNsigmaToVertex(0),
   fCutSigmaToVertexRequired(0),
+  fCutMaxDCAToVertexXY(0),
+  fCutMaxDCAToVertexZ(0),
+  fCutMinDCAToVertexXY(0),
+  fCutMinDCAToVertexZ(0),
+  fCutDCAToVertex2D(0),
   fPMin(0),
   fPMax(0),
   fPtMin(0),
@@ -190,16 +222,18 @@ AliESDtrackCuts::~AliESDtrackCuts()
     if (fhDXY[i])
       delete fhDXY[i];                     
     if (fhDZ[i])
-      delete fhDZ[i];                      
+      delete fhDZ[i];
+    if (fhDXYDZ[i])
+      delete fhDXYDZ[i];
     if (fhDXYvsDZ[i])
-      delete fhDXYvsDZ[i];                 
-    
+      delete fhDXYvsDZ[i];
+
     if (fhDXYNormalized[i])
       delete fhDXYNormalized[i];           
     if (fhDZNormalized[i])
       delete fhDZNormalized[i];
     if (fhDXYvsDZNormalized[i])
-      delete fhDXYvsDZNormalized[i];       
+      delete fhDXYvsDZNormalized[i];
     if (fhNSigmaToVertex[i])
       delete fhNSigmaToVertex[i];
     if (fhPt[i])
@@ -228,6 +262,9 @@ void AliESDtrackCuts::Init()
 
   fCutMaxChi2PerClusterTPC = 0;
   fCutMaxChi2PerClusterITS = 0;
+  
+  for (Int_t i = 0; i < 3; i++)
+       fCutClusterRequirementITS[i] = kOff;
 
   fCutMaxC11 = 0;
   fCutMaxC22 = 0;
@@ -238,10 +275,17 @@ void AliESDtrackCuts::Init()
   fCutAcceptKinkDaughters = 0;
   fCutRequireTPCRefit = 0;
   fCutRequireITSRefit = 0;
+  fCutRequireITSStandAlone = 0;
 
   fCutNsigmaToVertex = 0;
   fCutSigmaToVertexRequired = 0;
+  fCutMaxDCAToVertexXY = 0;
+  fCutMaxDCAToVertexZ = 0;
+  fCutDCAToVertex2D = 0;
+  fCutMinDCAToVertexXY = 0;
+  fCutMinDCAToVertexZ = 0;
 
+  
   fPMin = 0;
   fPMax = 0;
   fPtMin = 0;
@@ -275,6 +319,7 @@ void AliESDtrackCuts::Init()
 
     fhDXY[i] = 0;
     fhDZ[i] = 0;
+    fhDXYDZ[i] = 0;
     fhDXYvsDZ[i] = 0;
 
     fhDXYNormalized[i] = 0;
@@ -319,6 +364,9 @@ void AliESDtrackCuts::Copy(TObject &c) const
   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
 
+  for (Int_t i = 0; i < 3; i++)
+    target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
+
   target.fCutMaxC11 = fCutMaxC11;
   target.fCutMaxC22 = fCutMaxC22;
   target.fCutMaxC33 = fCutMaxC33;
@@ -328,9 +376,15 @@ void AliESDtrackCuts::Copy(TObject &c) const
   target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
   target.fCutRequireTPCRefit = fCutRequireTPCRefit;
   target.fCutRequireITSRefit = fCutRequireITSRefit;
+  target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
 
   target.fCutNsigmaToVertex = fCutNsigmaToVertex;
   target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
+  target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
+  target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
+  target.fCutDCAToVertex2D = fCutDCAToVertex2D;
+  target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
+  target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
 
   target.fPMin = fPMin;
   target.fPMax = fPMax;
@@ -365,6 +419,7 @@ void AliESDtrackCuts::Copy(TObject &c) const
 
     if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
     if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
+    if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
     if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
 
     if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
@@ -396,7 +451,6 @@ Long64_t AliESDtrackCuts::Merge(TCollection* list) {
   TIterator* iter = list->MakeIterator();
   TObject* obj;
 
-
   // collection of measured and generated histograms
   Int_t count = 0;
   while ((obj = iter->Next())) {
@@ -407,7 +461,7 @@ Long64_t AliESDtrackCuts::Merge(TCollection* list) {
 
     if (!entry->fHistogramsOn)
       continue;
-    
+
     for (Int_t i=0; i<2; i++) {
       
       fhNClustersITS[i]      ->Add(entry->fhNClustersITS[i]     );      
@@ -424,11 +478,12 @@ Long64_t AliESDtrackCuts::Merge(TCollection* list) {
                                                                    
       fhDXY[i]               ->Add(entry->fhDXY[i]              ); 
       fhDZ[i]                ->Add(entry->fhDZ[i]               ); 
-      fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          ); 
-                                                                   
-      fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    ); 
+      fhDXYDZ[i]             ->Add(entry->fhDXYDZ[i]          );
+      fhDXYvsDZ[i]           ->Add(entry->fhDXYvsDZ[i]          );
+
+      fhDXYNormalized[i]     ->Add(entry->fhDXYNormalized[i]    );
       fhDZNormalized[i]      ->Add(entry->fhDZNormalized[i]     );
-      fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]); 
+      fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
       fhNSigmaToVertex[i]    ->Add(entry->fhNSigmaToVertex[i]); 
 
       fhPt[i]                ->Add(entry->fhPt[i]); 
@@ -455,7 +510,7 @@ Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
   esdTrack->GetImpactParameters(b,bCov);
   
   if (bCov[0]<=0 || bCov[2]<=0) {
-    AliDebug(1, "Estimated b resolution lower or equal zero!");
+    AliDebugClass(1, "Estimated b resolution lower or equal zero!");
     bCov[0]=0; bCov[2]=0;
   }
   bRes[0] = TMath::Sqrt(bCov[0]);
@@ -469,7 +524,7 @@ Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
   //
-  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
+  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
   // Can this be expressed in a different way?
 
   if (bRes[0] == 0 || bRes[1] ==0)
@@ -477,13 +532,14 @@ Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
 
   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 
-  // stupid rounding problem screws up everything:
+  // work around precision problem
   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
-  if (TMath::Exp(-d * d / 2) < 1e-10)
+  // 1e-15 corresponds to nsigma ~ 7.7
+  if (TMath::Exp(-d * d / 2) < 1e-15)
     return 1000;
 
-  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
-  return d;
+  Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
+  return nSigma;
 }
 
 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
@@ -507,8 +563,8 @@ void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 }
 
 //____________________________________________________________________
-Bool_t
-AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
+Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) 
+{
   // 
   // figure out if the tracks survives all the track cuts defined
   //
@@ -533,16 +589,12 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   // fTracks.fP   //GetMass
   // fTracks.fKinkIndexes
 
-
   UInt_t status = esdTrack->GetStatus();
 
-  // dummy array
-  Int_t  fIdxInt[200];
-
   // getting quality parameters from the ESD track
-  Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
-  Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
-  
+  Int_t nClustersITS = esdTrack->GetITSclusters(0);
+  Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
+
   Float_t chi2PerClusterITS = -1;
   Float_t chi2PerClusterTPC = -1;
   if (nClustersITS!=0)
@@ -554,7 +606,27 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
 
   // getting the track to vertex parameters
   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
+      
+  Float_t b[2];
+  Float_t bCov[3];
+  esdTrack->GetImpactParameters(b,bCov);
+  if (bCov[0]<=0 || bCov[2]<=0) {
+    AliDebug(1, "Estimated b resolution lower or equal zero!");
+    bCov[0]=0; bCov[2]=0;
+  }
 
+  Float_t dcaToVertexXY = b[0];
+  Float_t dcaToVertexZ = b[1];
+
+  Float_t dcaToVertex = -1;
+  
+  if (fCutDCAToVertex2D)
+  {
+    dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
+  }
+  else
+    dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
+    
   // getting the kinematic variables of the track
   // (assuming the mass is known)
   Double_t p[3];
@@ -621,15 +693,35 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     cuts[17] = kTRUE;
   if((p[2] < fPzMin) || (p[2] > fPzMax))
     cuts[18] = kTRUE;
-  if((eta < fEtaMin) || (eta > fEtaMax)) 
+  if((eta < fEtaMin) || (eta > fEtaMax))
     cuts[19] = kTRUE;
   if((y < fRapMin) || (y > fRapMax)) 
     cuts[20] = kTRUE;
-
+  if (fCutDCAToVertex2D && dcaToVertex > 1)
+    cuts[21] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
+    cuts[22] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
+    cuts[23] = kTRUE;
+  if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
+    cuts[24] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
+    cuts[25] = kTRUE;
+  if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
+    cuts[26] = kTRUE;
+  
+  for (Int_t i = 0; i < 3; i++)
+    cuts[27+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
+  
+  if (fCutRequireITSStandAlone && ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)))
+    cuts[30]=kTRUE;
+  
   Bool_t cut=kFALSE;
   for (Int_t i=0; i<kNCuts; i++) 
-    if (cuts[i]) cut = kTRUE;
-  
+    if (cuts[i]) {cut = kTRUE;}
+
+
+
   //########################################################################
   // filling histograms
   if (fHistogramsOn) {
@@ -639,141 +731,133 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     
     for (Int_t i=0; i<kNCuts; i++) {
       if (cuts[i])
-       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
-      
+        fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
+
       for (Int_t j=i; j<kNCuts; j++) {
-       if (cuts[i] && cuts[j]) {
-         Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
-         Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
-         fhCutCorrelation->Fill(xC, yC);
-       }
+        if (cuts[i] && cuts[j]) {
+          Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
+          Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
+          fhCutCorrelation->Fill(xC, yC);
+        }
       }
     }
-    
-    fhNClustersITS[0]->Fill(nClustersITS);
-    fhNClustersTPC[0]->Fill(nClustersTPC);
-    fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
-    fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
-
-    fhC11[0]->Fill(extCov[0]);
-    fhC22[0]->Fill(extCov[2]);
-    fhC33[0]->Fill(extCov[5]);
-    fhC44[0]->Fill(extCov[9]);
-    fhC55[0]->Fill(extCov[14]);
-    
-    fhPt[0]->Fill(pt);
-    fhEta[0]->Fill(eta);
-
-    Float_t b[2];
-    Float_t bRes[2];
-    Float_t bCov[3];
-    esdTrack->GetImpactParameters(b,bCov);
-
-    if (bCov[0]<=0 || bCov[2]<=0) {
-      AliDebug(1, "Estimated b resolution lower or equal zero!");
-      bCov[0]=0; bCov[2]=0;
-    }
-    bRes[0] = TMath::Sqrt(bCov[0]);
-    bRes[1] = TMath::Sqrt(bCov[2]);
-
-    fhDZ[0]->Fill(b[1]);
-    fhDXY[0]->Fill(b[0]);
-    fhDXYvsDZ[0]->Fill(b[1],b[0]);
-
-    if (bRes[0]!=0 && bRes[1]!=0) {
-      fhDZNormalized[0]->Fill(b[1]/bRes[1]);
-      fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
-      fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
-      fhNSigmaToVertex[0]->Fill(nSigmaToVertex);
-    }
   }
 
-  //########################################################################
-  // cut the track!
-  if (cut) return kFALSE;
-
-  //########################################################################
-  // filling histograms after cut
-  if (fHistogramsOn) {
-    fhNClustersITS[1]->Fill(nClustersITS);
-    fhNClustersTPC[1]->Fill(nClustersTPC);
-    fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
-    fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
-
-    fhC11[1]->Fill(extCov[0]);
-    fhC22[1]->Fill(extCov[2]);
-    fhC33[1]->Fill(extCov[5]);
-    fhC44[1]->Fill(extCov[9]);
-    fhC55[1]->Fill(extCov[14]);
-
-    fhPt[1]->Fill(pt);
-    fhEta[1]->Fill(eta);
-    
-    Float_t b[2];
-    Float_t bRes[2];
-    Float_t bCov[3];
-    esdTrack->GetImpactParameters(b,bCov);
-    if (bCov[0]<=0 || bCov[2]<=0) {
-      AliDebug(1, "Estimated b resolution lower or equal zero!");
-      bCov[0]=0; bCov[2]=0;
-    }
-    bRes[0] = TMath::Sqrt(bCov[0]);
-    bRes[1] = TMath::Sqrt(bCov[2]);
+  // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
+  // the code is not in a function due to too many local variables that would need to be passed
 
-    fhDZ[1]->Fill(b[1]);
-    fhDXY[1]->Fill(b[0]);
-    fhDXYvsDZ[1]->Fill(b[1],b[0]);
+  for (Int_t id = 0; id < 2; id++)
+  {
+    // id = 0 --> before cut
+    // id = 1 --> after cut
 
-    if (bRes[0]!=0 && bRes[1]!=0)
+    if (fHistogramsOn)
     {
-      fhDZNormalized[1]->Fill(b[1]/bRes[1]);
-      fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
-      fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
-      fhNSigmaToVertex[1]->Fill(nSigmaToVertex);
+      fhNClustersITS[id]->Fill(nClustersITS);
+      fhNClustersTPC[id]->Fill(nClustersTPC);
+      fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
+      fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
+
+      fhC11[id]->Fill(extCov[0]);
+      fhC22[id]->Fill(extCov[2]);
+      fhC33[id]->Fill(extCov[5]);
+      fhC44[id]->Fill(extCov[9]);
+      fhC55[id]->Fill(extCov[14]);
+
+      fhPt[id]->Fill(pt);
+      fhEta[id]->Fill(eta);
+
+      Float_t bRes[2];
+      bRes[0] = TMath::Sqrt(bCov[0]);
+      bRes[1] = TMath::Sqrt(bCov[2]);
+
+      fhDZ[id]->Fill(b[1]);
+      fhDXY[id]->Fill(b[0]);
+      fhDXYDZ[id]->Fill(dcaToVertex);
+      fhDXYvsDZ[id]->Fill(b[1],b[0]);
+
+      if (bRes[0]!=0 && bRes[1]!=0) {
+        fhDZNormalized[id]->Fill(b[1]/bRes[1]);
+        fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
+        fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
+        fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
+      }
     }
+
+    // cut the track
+    if (cut)
+      return kFALSE;
   }
 
   return kTRUE;
 }
 
 //____________________________________________________________________
-TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
+Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
 {
-  //
-  // returns an array of all tracks that pass the cuts
-  //
-
-  TObjArray* acceptedTracks = new TObjArray();
-
-
-  // loop over esd tracks
-  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
-    AliESDtrack* track = esd->GetTrack(iTrack);
-    if (AcceptTrack(track))
-       acceptedTracks->Add(track);
+  // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
+  
+  switch (req)
+  {
+       case kOff:        return kTRUE;
+       case kNone:       return !clusterL1 && !clusterL2;
+       case kAny:        return clusterL1 || clusterL2;
+       case kFirst:      return clusterL1;
+       case kOnlyFirst:  return clusterL1 && !clusterL2;
+       case kSecond:     return clusterL2;
+       case kOnlySecond: return clusterL2 && !clusterL1;
+       case kBoth:       return clusterL1 && clusterL2;
   }
-  return acceptedTracks;
+  
+  return kFALSE;
 }
 
-
 //____________________________________________________________________
-Int_t AliESDtrackCuts::CountAcceptedTracks(AliESD* esd)
+AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
 {
+  // creates a TPC only track from the given esd track
+  // the track has to be deleted by the user
   //
-  // returns an the number of tracks that pass the cuts
+  // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
+  // there are only missing propagations here that are needed for old data
+  // this function will therefore become obsolete
   //
+  // adapted from code provided by CKB
 
-  Int_t count = 0;
+  if (!esd->GetPrimaryVertexTPC())
+    return 0; // No TPC vertex no TPC tracks
 
-  // loop over esd tracks
-  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
-    AliESDtrack* track = esd->GetTrack(iTrack);
+  if(!esd->GetPrimaryVertexTPC()->GetStatus())
+    return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
 
-    if (AcceptTrack(track))
-      count++;
+  AliESDtrack* track = esd->GetTrack(iTrack);
+  if (!track)
+    return 0;
+
+  AliESDtrack *tpcTrack = new AliESDtrack();
+
+  // This should have been done during the reconstruction
+  // fixed by Juri in r26675
+  // but recalculate for older data CKB
+  Float_t p[2],cov[3];
+  track->GetImpactParametersTPC(p,cov);
+  if(p[0]==0&&p[1]==0)
+    track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
+  // BKC
+
+  // only true if we have a tpc track
+  if (!track->FillTPCOnlyTrack(*tpcTrack))
+  {
+    delete tpcTrack;
+    return 0;
   }
 
-  return count;
+  // propagate to Vertex
+  // not needed for normal reconstructed ESDs...
+  // Double_t pTPC[2],covTPC[3];
+  // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
+
+  return tpcTrack;
 }
 
 //____________________________________________________________________
@@ -788,39 +872,26 @@ TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
 
   // loop over esd tracks
   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
-    AliESDtrack* track = esd->GetTrack(iTrack);
-
     if(bTPC){
       if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
+      if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
 
-      AliESDtrack *tpcTrack = new AliESDtrack();
-      bool bAdd = false;
-      Double_t pTPC[2],covTPC[3];
-      // This should have been done during the reconstruction
-      // fixed by Juri in r26675
-      // but recalculate for older data CKB 
-      Float_t p[2],cov[3];
-      track->GetImpactParametersTPC(p,cov);
-      if(p[0]==0&&p[1]==0){
-       track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
-      }
-      // BKC
-
-      if(track->FillTPCOnlyTrack(*tpcTrack)){ // only true if we have a tpc track
-       // propagate to Vertex
-       // not needed for normal reconstructed ESDs...
-       //      if(tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC))
-       if(AcceptTrack(tpcTrack)){
-         acceptedTracks->Add(tpcTrack);
-         bAdd = true;
-       }      
+      AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
+      if (!tpcTrack)
+        continue;
+
+      if (AcceptTrack(tpcTrack)) {
+        acceptedTracks->Add(tpcTrack);
       }
-      if(!bAdd)delete tpcTrack;
+      else
+        delete tpcTrack;
+    }
+    else
+    {
+      AliESDtrack* track = esd->GetTrack(iTrack);
+      if(AcceptTrack(track))
+        acceptedTracks->Add(track);
     }
-    else if(AcceptTrack(track)){// we cut by passing the original track
-      // default case
-      acceptedTracks->Add(track);
-    } 
   } 
   if(bTPC)acceptedTracks->SetOwner(kTRUE);
   return acceptedTracks;
@@ -870,41 +941,38 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
-   } 
+   }
 
   fhCutStatistics  ->SetLineColor(color);
   fhCutCorrelation ->SetLineColor(color);
   fhCutStatistics  ->SetLineWidth(2);
   fhCutCorrelation ->SetLineWidth(2);
 
-  Char_t str[256];
   for (Int_t i=0; i<2; i++) {
-    if (i==0) sprintf(str," ");
-    else sprintf(str,"_cut");
-
-    fhNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str)     ,"",8,-0.5,7.5);
-    fhNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str)     ,"",165,-0.5,164.5);
-    fhChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
-    fhChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
-
-    fhC11[i]                 = new  TH1F(Form("covMatrixDiagonal11%s",str),"",2000,0,20);
-    fhC22[i]                 = new  TH1F(Form("covMatrixDiagonal22%s",str),"",2000,0,20);
-    fhC33[i]                 = new  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,1);
-    fhC44[i]                 = new  TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
-    fhC55[i]                 = new  TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
-
-    fhDXY[i]                 = new  TH1F(Form("dXY%s",str)    ,"",500,-10,10);
-    fhDZ[i]                  = new  TH1F(Form("dZ%s",str)     ,"",500,-10,10);
-    fhDXYvsDZ[i]             = new  TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
-
-    fhDXYNormalized[i]       = new  TH1F(Form("dXYNormalized%s",str)    ,"",500,-10,10);
-    fhDZNormalized[i]        = new  TH1F(Form("dZNormalized%s",str)     ,"",500,-10,10);
-    fhDXYvsDZNormalized[i]   = new  TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
-
-    fhNSigmaToVertex[i]      = new  TH1F(Form("nSigmaToVertex%s",str),"",500,0,50);
-
-    fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
-    fhEta[i]                 = new TH1F(Form("eta%s",str)     ,"#eta distribution;#eta",40,-2.0,2.0);
+    fhNClustersITS[i]        = new TH1F("nClustersITS"    ,"",8,-0.5,7.5);
+    fhNClustersTPC[i]        = new TH1F("nClustersTPC"     ,"",165,-0.5,164.5);
+    fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
+    fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
+
+    fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
+    fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
+    fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
+    fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
+    fhC55[i]                 = new TH1F("covMatrixDiagonal55","",1000,0,5);
+
+    fhDXY[i]                 = new TH1F("dXY"    ,"",500,-10,10);
+    fhDZ[i]                  = new TH1F("dZ"     ,"",500,-10,10);
+    fhDXYDZ[i]               = new TH1F("dXYDZ"  ,"",500,0,10);
+    fhDXYvsDZ[i]             = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
+
+    fhDXYNormalized[i]       = new TH1F("dXYNormalized"    ,"",500,-10,10);
+    fhDZNormalized[i]        = new TH1F("dZNormalized"     ,"",500,-10,10);
+    fhDXYvsDZNormalized[i]   = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
+
+    fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
+
+    fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
+    fhEta[i]                 = new TH1F("eta"     ,"#eta distribution;#eta",40,-2.0,2.0);
     
     fhNClustersITS[i]->SetTitle("n ITS clusters");
     fhNClustersTPC[i]->SetTitle("n TPC clusters");
@@ -917,15 +985,16 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
     fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
     fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
 
-    fhDXY[i]->SetTitle("transverse impact parameter");
-    fhDZ[i]->SetTitle("longitudinal impact parameter");
-    fhDXYvsDZ[i]->SetTitle("longitudinal impact parameter");
-    fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
+    fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
+    fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
+    fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2)  (cm)");
+    fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
+    fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
 
-    fhDXYNormalized[i]->SetTitle("normalized trans impact par");
-    fhDZNormalized[i]->SetTitle("normalized long impact par");
-    fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par");
-    fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
+    fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
+    fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
+    fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
+    fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
     fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
 
     fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
@@ -940,7 +1009,8 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
     fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
 
     fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
-    fhDZ[i]->SetLineColor(color);   fhDZ[i]->SetLineWidth(2);
+    fhDZ[i]->SetLineColor(color);    fhDZ[i]->SetLineWidth(2);
+    fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
 
     fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
     fhDZNormalized[i]->SetLineColor(color);    fhDZNormalized[i]->SetLineWidth(2);
@@ -950,7 +1020,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
   // The number of sigmas to the vertex is per definition gaussian
   ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
   ffDTheoretical->SetParameter(0,1);
-  
+
   TH1::AddDirectory(oldStatus);
 }
 
@@ -973,41 +1043,37 @@ Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
 
-  Char_t str[5];
   for (Int_t i=0; i<2; i++) {
     if (i==0)
     {
       gDirectory->cd("before_cuts");
-      str[0] = 0;
     }
     else
-    {
       gDirectory->cd("after_cuts");
-      sprintf(str,"_cut");
-    }
 
-    fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("nClustersITS%s",str)     ));
-    fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("nClustersTPC%s",str)     ));
-    fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("chi2PerClusterITS%s",str)));
-    fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("chi2PerClusterTPC%s",str)));
+    fhNClustersITS[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS"     ));
+    fhNClustersTPC[i]      = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC"     ));
+    fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
+    fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
 
-    fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("covMatrixDiagonal11%s",str)));
-    fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("covMatrixDiagonal22%s",str)));
-    fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("covMatrixDiagonal33%s",str)));
-    fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("covMatrixDiagonal44%s",str)));
-    fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("covMatrixDiagonal55%s",str)));
+    fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
+    fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
+    fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
+    fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
+    fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
 
-    fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get(Form("dXY%s",str)    ));
-    fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get(Form("dZ%s",str)     ));
-    fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get(Form("dXYvsDZ%s",str)));
+    fhDXY[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXY"    ));
+    fhDZ[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZ"     ));
+    fhDXYDZ[i] =   dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
+    fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
 
-    fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get(Form("dXYNormalized%s",str)    ));
-    fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get(Form("dZNormalized%s",str)     ));
-    fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get(Form("dXYvsDZNormalized%s",str)));
-    fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("nSigmaToVertex%s",str)));
+    fhDXYNormalized[i] =     dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized"    ));
+    fhDZNormalized[i] =      dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized"     ));
+    fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
+    fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
 
-    fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
-    fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("eta%s",str)));
+    fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
+    fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
 
     gDirectory->cd("../");
   }
@@ -1063,6 +1129,7 @@ void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
 
     fhDXY[i]                 ->Write();
     fhDZ[i]                  ->Write();
+    fhDXYDZ[i]               ->Write();
     fhDXYvsDZ[i]             ->Write();
 
     fhDXYNormalized[i]       ->Write();
@@ -1208,3 +1275,38 @@ void AliESDtrackCuts::DrawHistograms()
   fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
 }
 
+Float_t AliESDtrackCuts::GetMinNsigmaToVertex() const
+{
+  // deprecated, please use GetMaxNsigmaToVertex
+
+  Printf("WARNING: AliESDtrackCuts::GetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use GetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
+
+  return GetMaxNsigmaToVertex();
+}
+
+void AliESDtrackCuts::SetMinNsigmaToVertex(Float_t sigma)
+{
+  // deprecated, will be removed in next release
+
+  Printf("WARNING: AliESDtrackCuts::SetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use SetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
+  
+  SetMaxNsigmaToVertex(sigma);
+}
+
+void AliESDtrackCuts::SetAcceptKingDaughters(Bool_t b)
+{
+  // deprecated, will be removed in next release
+
+  Printf("WARNING: AliESDtrackCuts::SetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use SetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
+  
+  SetAcceptKinkDaughters(b);
+}
+
+Bool_t AliESDtrackCuts::GetAcceptKingDaughters() const
+{
+  // deprecated, will be removed in next release
+
+  Printf("WARNING: AliESDtrackCuts::GetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use GetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
+  
+  return GetAcceptKinkDaughters();
+}