added CountAcceptedTracks function
[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx
index e71833ce7181af4aa1e05273a688f23272fb8bdb..98ead8e1349c800487d0af42dabdb4851fbb2958 100644 (file)
@@ -1,5 +1,10 @@
 #include "AliESDtrackCuts.h"
 
+
+#include <AliESDtrack.h>
+#include <AliESD.h>
+#include <AliLog.h>
+
 //____________________________________________________________________
 ClassImp(AliESDtrackCuts)
 
@@ -246,6 +251,8 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   // 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
@@ -255,6 +262,8 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
   Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
   
+
+
   Float_t chi2PerClusterITS = -1;
   Float_t chi2PerClusterTPC = -1;
   if (nClustersITS!=0)
@@ -271,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)
@@ -293,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.;
@@ -441,11 +464,11 @@ AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
   //
 
   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);
   }
@@ -453,6 +476,27 @@ AliESDtrackCuts::GetAcceptedTracks(AliESD* esd)
   return acceptedTracks;
 }
 
+//____________________________________________________________________
+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);
+
+    if (AcceptTrack(track))
+      count++;
+  }
+
+  return count;
+}
+
 //____________________________________________________________________
  void AliESDtrackCuts::DefineHistograms(Int_t color) {
    //