added CountAcceptedTracks function
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
index 6b0b1051a17275c78fae85084b46f5521ae16455..98ead8e1349c800487d0af42dabdb4851fbb2958 100644 (file)
@@ -1,18 +1,52 @@
 #include "AliESDtrackCuts.h"
 
-#include <Riostream.h>
+
+#include <AliESDtrack.h>
+#include <AliESD.h>
+#include <AliLog.h>
 
 //____________________________________________________________________
-ClassImp(AliESDtrackCuts);
+ClassImp(AliESDtrackCuts)
+
+// Cut names
+const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
+ "require TPC refit",
+ "require ITS refit",
+ "n clusters TPC",
+ "n clusters ITS",
+ "#Chi^{2}/clusters TPC",
+ "#Chi^{2}/clusters ITS",
+ "cov 11",
+ "cov 22",
+ "cov 33",
+ "cov 44",
+ "cov 55",
+ "trk-to-vtx",
+ "trk-to-vtx failed",
+ "kink daughters",
+
+ "p",
+ "p_{T}",
+ "p_{x}",
+ "p_{y}",
+ "p_{z}",
+ "y",
+ "eta"
+};
 
 //____________________________________________________________________
-AliESDtrackCuts::AliESDtrackCuts() {
+AliESDtrackCuts::AliESDtrackCuts()
+{
+  //
+  // constructor
+  //
+
+  Init();
 
   //##############################################################################
   // setting default cuts
-
   SetMinNClustersTPC();
-  SetMinNClustersITS();            
+  SetMinNClustersITS();
   SetMaxChi2PerClusterTPC();
   SetMaxChi2PerClusterITS();                               
   SetMaxCovDiagonalElements();                                     
@@ -30,51 +64,182 @@ AliESDtrackCuts::AliESDtrackCuts() {
   SetRapRange();
 
   SetHistogramsOn();
+}
 
-  // set the cut names
-  fCutNames[0]  = "require TPC refit";
-  fCutNames[1]  = "require ITS refit";  
-  fCutNames[2]  = "n clusters TPC";
-  fCutNames[3]  = "n clusters ITS";
-  fCutNames[4]  = "#Chi^{2}/clusters TPC";
-  fCutNames[5]  = "#Chi^{2}/clusters ITS";
-  fCutNames[6]  = "cov 11";
-  fCutNames[7]  = "cov 22";
-  fCutNames[8]  = "cov 33";
-  fCutNames[9]  = "cov 44";
-  fCutNames[10] = "cov 55";
-  fCutNames[11] = "trk-to-vtx";
-  fCutNames[12] = "trk-to-vtx failed";
-  fCutNames[13] = "kink daughters";
-
-  fCutNames[14] = "p";
-  fCutNames[15] = "p_{T}";
-  fCutNames[16] = "p_{x}";
-  fCutNames[17] = "p_{y}";
-  fCutNames[18] = "p_{z}";
-  fCutNames[19] = "y";
-  fCutNames[20] = "eta";
+//_____________________________________________________________________________
+AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : TObject(c)
+{
+  //
+  // copy constructor
+  //
 
+  ((AliESDtrackCuts &) c).Copy(*this);
 }
 
-//____________________________________________________________________
-Bool_t 
-AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {
+AliESDtrackCuts::~AliESDtrackCuts()
+{
   //
-  // re-calculate the track-to-vertex
-  esdTrack->RelateToVertex(esdVtx, field, 1e99);
-  
-  return AcceptTrack(esdTrack);
+  // destructor
+  //
+
+  // ## TODO to be implemented
 }
 
-//____________________________________________________________________
-Bool_t 
-AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) {
-  
-  AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex");
-  esdTrack->RelateToVertex(esdVtx, field, 1e99);
-  return AcceptTrack(esdTrack);
-} 
+void AliESDtrackCuts::Init()
+{
+  //
+  // sets everything to zero
+  //
+
+  fCutMinNClusterTPC = 0;
+  fCutMinNClusterITS = 0;
+
+  fCutMaxChi2PerClusterTPC = 0;
+  fCutMaxChi2PerClusterITS = 0;
+
+  fCutMaxC11 = 0;
+  fCutMaxC22 = 0;
+  fCutMaxC33 = 0;
+  fCutMaxC44 = 0;
+  fCutMaxC55 = 0;
+
+  fCutAcceptKinkDaughters = 0;
+  fCutRequireTPCRefit = 0;
+  fCutRequireITSRefit = 0;
+
+  fCutNsigmaToVertex = 0;
+  fCutSigmaToVertexRequired = 0;
+
+  fPMin = 0;
+  fPMax = 0;
+  fPtMin = 0;
+  fPtMax = 0;
+  fPxMin = 0;
+  fPxMax = 0;
+  fPyMin = 0;
+  fPyMax = 0;
+  fPzMin = 0;
+  fPzMax = 0;
+  fEtaMin = 0;
+  fEtaMax = 0;
+  fRapMin = 0;
+  fRapMax = 0;
+
+  fHistogramsOn = kFALSE;
+
+  for (Int_t i=0; i<2; ++i)
+  {
+    fhNClustersITS[i] = 0;
+    fhNClustersTPC[i] = 0;
+
+    fhChi2PerClusterITS[i] = 0;
+    fhChi2PerClusterTPC[i] = 0;
+
+    fhC11[i] = 0;
+    fhC22[i] = 0;
+    fhC33[i] = 0;
+    fhC44[i] = 0;
+    fhC55[i] = 0;
+
+    fhDXY[i] = 0;
+    fhDZ[i] = 0;
+    fhDXYvsDZ[i] = 0;
+
+    fhDXYNormalized[i] = 0;
+    fhDZNormalized[i] = 0;
+    fhDXYvsDZNormalized[i] = 0;
+  }
+
+  fhCutStatistics = 0;
+  fhCutCorrelation = 0;
+}
+
+//_____________________________________________________________________________
+AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+void AliESDtrackCuts::Copy(TObject &c) const
+{
+  //
+  // Copy function
+  //
+
+  AliESDtrackCuts& target = (AliESDtrackCuts &) c;
+
+  target.Init();
+
+  target.fCutMinNClusterTPC = fCutMinNClusterTPC;
+  target.fCutMinNClusterITS = fCutMinNClusterITS;
+
+  target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
+  target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
+
+  target.fCutMaxC11 = fCutMaxC11;
+  target.fCutMaxC22 = fCutMaxC22;
+  target.fCutMaxC33 = fCutMaxC33;
+  target.fCutMaxC44 = fCutMaxC44;
+  target.fCutMaxC55 = fCutMaxC55;
+
+  target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
+  target.fCutRequireTPCRefit = fCutRequireTPCRefit;
+  target.fCutRequireITSRefit = fCutRequireITSRefit;
+
+  target.fCutNsigmaToVertex = fCutNsigmaToVertex;
+  target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
+
+  target.fPMin = fPMin;
+  target.fPMax = fPMax;
+  target.fPtMin = fPtMin;
+  target.fPtMax = fPtMax;
+  target.fPxMin = fPxMin;
+  target.fPxMax = fPxMax;
+  target.fPyMin = fPyMin;
+  target.fPyMax = fPyMax;
+  target.fPzMin = fPzMin;
+  target.fPzMax = fPzMax;
+  target.fEtaMin = fEtaMin;
+  target.fEtaMax = fEtaMax;
+  target.fRapMin = fRapMin;
+  target.fRapMax = fRapMax;
+
+  target.fHistogramsOn = fHistogramsOn;
+
+  for (Int_t i=0; i<2; ++i)
+  {
+    if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
+    if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
+
+    if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
+    if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
+
+    if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
+    if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
+    if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
+    if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
+    if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
+
+    if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
+    if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
+    if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
+
+    if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
+    if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
+    if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
+  }
+
+  if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
+  if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
+
+  TObject::Copy(c);
+}
 
 //____________________________________________________________________
 Bool_t 
@@ -82,13 +247,23 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   // 
   // figure out if the tracks survives all the track cuts defined
   //
+  // the different quality parameter and kinematic values are first
+  // retrieved from the track. then it is found out what cuts the
+  // track did not survive and finally the cuts are imposed.
+
+
 
   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);
   
+
+
   Float_t chi2PerClusterITS = -1;
   Float_t chi2PerClusterTPC = -1;
   if (nClustersITS!=0)
@@ -97,7 +272,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);  
 
   Double_t extCov[15];
-  esdTrack->GetExternalCovariance(extCov);  
+  esdTrack->GetExternalCovariance(extCov);
 
   // getting the track to vertex parameters
   Float_t b[2];
@@ -105,19 +280,32 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   Float_t bCov[3];
   esdTrack->GetImpactParameters(b,bCov);    
   if (bCov[0]<=0 || bCov[2]<=0) {
-    AliDebug(1, "Estimated b resolution zero!");
-    bCov[0]=0; bCov[1]=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]);
 
-  // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  // -----------------------------------
+  // How to get to a n-sigma cut?
   //
-  // this is not correct - it will not give n sigma!!!
-  // 
+  // The accumulated statistics from 0 to d is
+  //
+  // ->  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)
+  // Can this be expressed in a different way?
+  //
+  //
+  // FIX: I don't think this is correct!!! Keeping d as n_sigma for now...
+
   Float_t nSigmaToVertex = -1;
-  if (bRes[0]!=0 && bRes[1]!=0)
-    nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));  
+  if (bRes[0]!=0 && bRes[1]!=0) {
+    Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
+    nSigmaToVertex = d;//TMath::Sqrt(2)*(TMath::ErfInverse(1 - TMath::Exp(0.5*(-d*d))));
+    // JF solution: nSigmaToVertex = TMath::ErfInverse(TMath::Sqrt(2.0/TMath::Pi()) * TMath::Erf(d / TMath::Sqrt(2))) * TMath::Sqrt(2);
+  }
 
   // getting the kinematic variables of the track 
   // (assuming the mass is known)
@@ -127,6 +315,7 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
   Float_t energy   = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
 
+
   //y-eta related calculations
   Float_t eta = -100.;
   Float_t y   = -100.;
@@ -139,38 +328,38 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   //########################################################################
   // cut the track?
   
-  Bool_t cuts[fNCuts];
-  for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
+  Bool_t cuts[kNCuts];
+  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
   
   // track quality cuts
-  if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
+  if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
     cuts[0]=kTRUE;
-  if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
+  if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
     cuts[1]=kTRUE;
-  if (nClustersTPC<fCut_MinNClusterTPC) 
+  if (nClustersTPC<fCutMinNClusterTPC) 
     cuts[2]=kTRUE;
-  if (nClustersITS<fCut_MinNClusterITS) 
+  if (nClustersITS<fCutMinNClusterITS) 
     cuts[3]=kTRUE;
-  if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC) 
+  if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC) 
     cuts[4]=kTRUE; 
-  if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS) 
+  if (chi2PerClusterITS>fCutMaxChi2PerClusterITS) 
     cuts[5]=kTRUE;
-  if (extCov[0]  > fCut_MaxC11) 
+  if (extCov[0]  > fCutMaxC11) 
     cuts[6]=kTRUE;  
-  if (extCov[2]  > fCut_MaxC22) 
+  if (extCov[2]  > fCutMaxC22) 
     cuts[7]=kTRUE;  
-  if (extCov[5]  > fCut_MaxC33) 
+  if (extCov[5]  > fCutMaxC33) 
     cuts[8]=kTRUE;  
-  if (extCov[9]  > fCut_MaxC44) 
+  if (extCov[9]  > fCutMaxC44) 
     cuts[9]=kTRUE;  
-  if (extCov[14]  > fCut_MaxC55) 
+  if (extCov[14]  > fCutMaxC55) 
     cuts[10]=kTRUE;  
-  if (nSigmaToVertex > fCut_NsigmaToVertex) 
+  if (nSigmaToVertex > fCutNsigmaToVertex) 
     cuts[11] = kTRUE;
   // if n sigma could not be calculated
-  if (nSigmaToVertex<0 && fCut_SigmaToVertexRequired)   
+  if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)   
     cuts[12]=kTRUE;
-  if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
+  if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
     cuts[13]=kTRUE;
   // track kinematics cut
   if((momentum < fPMin) || (momentum > fPMax)) 
@@ -189,50 +378,50 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     cuts[20] = kTRUE;
 
   Bool_t cut=kFALSE;
-  for (Int_t i=0; i<fNCuts; i++) 
+  for (Int_t i=0; i<kNCuts; i++) 
     if (cuts[i]) cut = kTRUE;
   
   //########################################################################
   // filling histograms
   if (fHistogramsOn) {
-    hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
+    fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
     
     if (cut)
-      hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
+      fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
     
-    for (Int_t i=0; i<fNCuts; i++) {
+    for (Int_t i=0; i<kNCuts; i++) {
       if (cuts[i])
-       hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
+       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
       
-      for (Int_t j=i; j<fNCuts; j++) {
+      for (Int_t j=i; j<kNCuts; j++) {
        if (cuts[i] && cuts[j]) {
-         Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
-         Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
-         hCutCorrelation->Fill(x,y);
+         Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
+         Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
+         fhCutCorrelation->Fill(x,y);
        }
       }
     }
     
 
-    hNClustersITS[0]->Fill(nClustersITS);        
-    hNClustersTPC[0]->Fill(nClustersTPC);        
-    hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
-    hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);   
+    fhNClustersITS[0]->Fill(nClustersITS);        
+    fhNClustersTPC[0]->Fill(nClustersTPC);        
+    fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
+    fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);   
     
-    hC11[0]->Fill(extCov[0]);                 
-    hC22[0]->Fill(extCov[2]);                 
-    hC33[0]->Fill(extCov[5]);                 
-    hC44[0]->Fill(extCov[9]);                                  
-    hC55[0]->Fill(extCov[14]);                                  
+    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]);                                  
     
-    hDZ[0]->Fill(b[1]);     
-    hDXY[0]->Fill(b[0]);    
-    hDXYvsDZ[0]->Fill(b[1],b[0]);
+    fhDZ[0]->Fill(b[1]);     
+    fhDXY[0]->Fill(b[0]);    
+    fhDXYvsDZ[0]->Fill(b[1],b[0]);
 
     if (bRes[0]!=0 && bRes[1]!=0) {
-      hDZNormalized[0]->Fill(b[1]/bRes[1]);     
-      hDXYNormalized[0]->Fill(b[0]/bRes[0]);    
-      hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[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]);
     }
   }
 
@@ -243,24 +432,24 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   //########################################################################  
   // filling histograms after cut
   if (fHistogramsOn) {
-    hNClustersITS[1]->Fill(nClustersITS);        
-    hNClustersTPC[1]->Fill(nClustersTPC);        
-    hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
-    hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);   
+    fhNClustersITS[1]->Fill(nClustersITS);        
+    fhNClustersTPC[1]->Fill(nClustersTPC);        
+    fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
+    fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);   
     
-    hC11[1]->Fill(extCov[0]);                 
-    hC22[1]->Fill(extCov[2]);                 
-    hC33[1]->Fill(extCov[5]);                 
-    hC44[1]->Fill(extCov[9]);                                  
-    hC55[1]->Fill(extCov[14]);                                  
+    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]);                                  
     
-    hDZ[1]->Fill(b[1]);     
-    hDXY[1]->Fill(b[0]);    
-    hDXYvsDZ[1]->Fill(b[1],b[0]);
+    fhDZ[1]->Fill(b[1]);     
+    fhDXY[1]->Fill(b[0]);    
+    fhDXYvsDZ[1]->Fill(b[1],b[0]);
 
-    hDZNormalized[1]->Fill(b[1]/bRes[1]);     
-    hDXYNormalized[1]->Fill(b[0]/bRes[0]);    
-    hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
+    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]);
   }
   
   return kTRUE;
@@ -268,146 +457,157 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
 
 //____________________________________________________________________
 TObjArray*
-AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) {
-  
+AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
+{
+  //
   // returns an array of all tracks that pass the cuts
-  fAcceptedTracks->Clear();
-  
+  //
+
+  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)) fAcceptedTracks->Add(track);
+
+    if (AcceptTrack(track))
+      acceptedTracks->Add(track);
   }
 
-  return fAcceptedTracks;
+  return acceptedTracks;
 }
 
 //____________________________________________________________________
-void 
-AliESDtrackCuts::DefineHistograms(Int_t color) {
+Int_t
+AliESDtrackCuts::CountAcceptedTracks(AliESD* esd)
+{
+  //
+  // returns an the number of tracks that pass the cuts
+  //
+
+  Int_t count = 0;
+
+  // loop over esd tracks
+  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
+    AliESDtrack* track = esd->GetTrack(iTrack);
 
-  fHistogramsOn=kTRUE;
+    if (AcceptTrack(track))
+      count++;
+  }
 
-  //###################################################################################
-  // defining histograms
+  return count;
+}
 
-  hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
+//____________________________________________________________________
+ void AliESDtrackCuts::DefineHistograms(Int_t color) {
+   // 
+   // diagnostics histograms are defined
+   // 
 
-  hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
-  hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
+   fHistogramsOn=kTRUE;
 
-  hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
-  
-  for (Int_t i=0; i<fNCuts; i++) {
-    hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
-    hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
-    hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
-  } 
+   //###################################################################################
+   // defining histograms
 
-  hCutStatistics  ->SetLineColor(color);
-  hCutCorrelation ->SetLineColor(color);
-  hCutStatistics  ->SetLineWidth(2);
-  hCutCorrelation ->SetLineWidth(2);
-
-
-  hNClustersITS        = new TH1F*[2];
-  hNClustersTPC        = new TH1F*[2];
-  hChi2PerClusterITS   = new TH1F*[2];
-  hChi2PerClusterTPC   = new TH1F*[2];
-                      
-  hC11                 = new TH1F*[2];
-  hC22                 = new TH1F*[2];
-  hC33                 = new TH1F*[2];
-  hC44                 = new TH1F*[2];
-  hC55                 = new TH1F*[2];
-  
-  hDXY                 = new TH1F*[2];
-  hDZ                 = new TH1F*[2];
-  hDXYvsDZ            = new TH2F*[2];
+   fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
+
+   fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
+   fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 
-  hDXYNormalized       = new TH1F*[2];
-  hDZNormalized        = new TH1F*[2];
-  hDXYvsDZNormalized   = new TH2F*[2];
+   fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
+  
+   for (Int_t i=0; i<kNCuts; i++) {
+     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");
 
-    hNClustersITS[i]        = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
-    hNClustersTPC[i]        = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
-    hChi2PerClusterITS[i]   = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
-    hChi2PerClusterTPC[i]   = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
-    
-    hC11[i]                 = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
-    hC22[i]                 = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
-    hC33[i]                 = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
-    hC44[i]                 = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
-    hC55[i]                 = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
-    
-    hDXY[i]                 = new TH1F(Form("dXY%s",str),"",500,-10,10);
-    hDZ[i]                  = new TH1F(Form("dZ%s",str),"",500,-10,10);
-    hDXYvsDZ[i]             = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
-
-    hDXYNormalized[i]       = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
-    hDZNormalized[i]        = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
-    hDXYvsDZNormalized[i]   = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
-
-
-    hNClustersITS[i]        ->SetXTitle("n ITS clusters");  
-    hNClustersTPC[i]        ->SetXTitle("n TPC clusters"); 
-    hChi2PerClusterITS[i]   ->SetXTitle("#Chi^{2} per ITS cluster"); 
-    hChi2PerClusterTPC[i]   ->SetXTitle("#Chi^{2} per TPC cluster"); 
-    
-    hC11[i]                 ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]"); 
-    hC22[i]                 ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]"); 
-    hC33[i]                 ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}"); 
-    hC44[i]                 ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}"); 
-    hC55[i]                 ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]"); 
-   
-    hDXY[i]                 ->SetXTitle("transverse impact parameter"); 
-    hDZ[i]                  ->SetXTitle("longitudinal impact parameter"); 
-    hDXYvsDZ[i]             ->SetXTitle("longitudinal impact parameter"); 
-    hDXYvsDZ[i]             ->SetYTitle("transverse impact parameter"); 
-
-    hDXYNormalized[i]       ->SetXTitle("normalized trans impact par"); 
-    hDZNormalized[i]        ->SetXTitle("normalized long impact par"); 
-    hDXYvsDZNormalized[i]   ->SetXTitle("normalized long impact par"); 
-    hDXYvsDZNormalized[i]   ->SetYTitle("normalized trans impact par"); 
-
-    hNClustersITS[i]        ->SetLineColor(color);   hNClustersITS[i]        ->SetLineWidth(2);
-    hNClustersTPC[i]        ->SetLineColor(color);   hNClustersTPC[i]        ->SetLineWidth(2);
-    hChi2PerClusterITS[i]   ->SetLineColor(color);   hChi2PerClusterITS[i]   ->SetLineWidth(2);
-    hChi2PerClusterTPC[i]   ->SetLineColor(color);   hChi2PerClusterTPC[i]   ->SetLineWidth(2);
-                                                                                             
-    hC11[i]                 ->SetLineColor(color);   hC11[i]                 ->SetLineWidth(2);
-    hC22[i]                 ->SetLineColor(color);   hC22[i]                 ->SetLineWidth(2);
-    hC33[i]                 ->SetLineColor(color);   hC33[i]                 ->SetLineWidth(2);
-    hC44[i]                 ->SetLineColor(color);   hC44[i]                 ->SetLineWidth(2);
-    hC55[i]                 ->SetLineColor(color);   hC55[i]                 ->SetLineWidth(2);
-                                                                                             
-    hDXY[i]                 ->SetLineColor(color);   hDXY[i]                 ->SetLineWidth(2);
-    hDZ[i]                  ->SetLineColor(color);   hDZ[i]                  ->SetLineWidth(2);
-                                                    
-    hDXYNormalized[i]       ->SetLineColor(color);   hDXYNormalized[i]       ->SetLineWidth(2);
-    hDZNormalized[i]        ->SetLineColor(color);   hDZNormalized[i]        ->SetLineWidth(2);
-
+    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),"",1000,0,5);
+    fhC22[i]                 = new  TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
+    fhC33[i]                 = new  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
+    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);
+
+
+    fhNClustersITS[i]->SetXTitle("n ITS clusters");
+    fhNClustersTPC[i]->SetXTitle("n TPC clusters");
+    fhChi2PerClusterITS[i]->SetXTitle("#Chi^{2} per ITS cluster");
+    fhChi2PerClusterTPC[i]->SetXTitle("#Chi^{2} per TPC cluster");
+
+    fhC11[i]->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
+    fhC22[i]->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
+    fhC33[i]->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
+    fhC44[i]->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
+    fhC55[i]->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
+
+    fhDXY[i]->SetXTitle("transverse impact parameter");
+    fhDZ[i]->SetXTitle("longitudinal impact parameter");
+    fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter");
+    fhDXYvsDZ[i]->SetYTitle("transverse impact parameter");
+
+    fhDXYNormalized[i]->SetXTitle("normalized trans impact par");
+    fhDZNormalized[i]->SetXTitle("normalized long impact par");
+    fhDXYvsDZNormalized[i]->SetXTitle("normalized long impact par"); 
+    fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par");
+
+    fhNClustersITS[i]->SetLineColor(color);   fhNClustersITS[i]->SetLineWidth(2);
+    fhNClustersTPC[i]->SetLineColor(color);   fhNClustersTPC[i]->SetLineWidth(2);
+    fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
+    fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
+
+    fhC11[i]->SetLineColor(color);   fhC11[i]->SetLineWidth(2);
+    fhC22[i]->SetLineColor(color);   fhC22[i]->SetLineWidth(2);
+    fhC33[i]->SetLineColor(color);   fhC33[i]->SetLineWidth(2);
+    fhC44[i]->SetLineColor(color);   fhC44[i]->SetLineWidth(2);
+    fhC55[i]->SetLineColor(color);   fhC55[i]->SetLineWidth(2);
+
+    fhDXY[i]->SetLineColor(color);   fhDXY[i]->SetLineWidth(2);
+    fhDZ[i]->SetLineColor(color);   fhDZ[i]->SetLineWidth(2);
+
+    fhDXYNormalized[i]->SetLineColor(color);   fhDXYNormalized[i]->SetLineWidth(2);
+    fhDZNormalized[i]->SetLineColor(color);   fhDZNormalized[i]->SetLineWidth(2);
   }
 }
 
 //____________________________________________________________________
 void 
-AliESDtrackCuts::Print() {
+AliESDtrackCuts::Print(const Option_t*) const {
+  //
+  // print method - still to be implemented
+  //
 
   AliInfo("AliESDtrackCuts...");
 }
 
 
 //____________________________________________________________________
-void 
-AliESDtrackCuts::SaveHistograms(Char_t* dir) {
+void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
+  // 
+  // saves the histograms in a directory (dir)
+  // 
+
   
   if (!fHistogramsOn) {
     AliDebug(0, "Histograms not on - cannot save histograms!!!");
@@ -420,8 +620,8 @@ AliESDtrackCuts::SaveHistograms(Char_t* dir) {
   gDirectory->mkdir("before_cuts");
   gDirectory->mkdir("after_cuts");
  
-  hCutStatistics->Write();
-  hCutCorrelation->Write();
+  fhCutStatistics->Write();
+  fhCutCorrelation->Write();
 
   for (Int_t i=0; i<2; i++) {
     if (i==0)
@@ -429,24 +629,24 @@ AliESDtrackCuts::SaveHistograms(Char_t* dir) {
     else
       gDirectory->cd("after_cuts");
     
-    hNClustersITS[i]        ->Write();
-    hNClustersTPC[i]        ->Write();
-    hChi2PerClusterITS[i]   ->Write();
-    hChi2PerClusterTPC[i]   ->Write();
+    fhNClustersITS[i]        ->Write();
+    fhNClustersTPC[i]        ->Write();
+    fhChi2PerClusterITS[i]   ->Write();
+    fhChi2PerClusterTPC[i]   ->Write();
     
-    hC11[i]                 ->Write();
-    hC22[i]                 ->Write();
-    hC33[i]                 ->Write();
-    hC44[i]                 ->Write();
-    hC55[i]                 ->Write();
-
-    hDXY[i]                 ->Write();
-    hDZ[i]                  ->Write();
-    hDXYvsDZ[i]             ->Write();
+    fhC11[i]                 ->Write();
+    fhC22[i]                 ->Write();
+    fhC33[i]                 ->Write();
+    fhC44[i]                 ->Write();
+    fhC55[i]                 ->Write();
+
+    fhDXY[i]                 ->Write();
+    fhDZ[i]                  ->Write();
+    fhDXYvsDZ[i]             ->Write();
     
-    hDXYNormalized[i]       ->Write();
-    hDZNormalized[i]        ->Write();
-    hDXYvsDZNormalized[i]   ->Write();
+    fhDXYNormalized[i]       ->Write();
+    fhDZNormalized[i]        ->Write();
+    fhDXYvsDZNormalized[i]   ->Write();
 
     gDirectory->cd("../");
   }