]> 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 13ced38a224e5f6bb99dd92140487ae634cbbe49..1a07eee202f996a256b3f4ad36c7ffbd478e43c6 100644 (file)
 #include "AliESDtrackCuts.h"
 
 #include <AliESDtrack.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)
@@ -34,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",
@@ -49,10 +52,18 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  "p_{x}",
  "p_{y}",
  "p_{z}",
- "y",
  "eta",
- "trk-to-vtx dca absolute",
- "trk-to-vtx dca xy absolute"
+ "y",
+ "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"
 };
 
 //____________________________________________________________________
@@ -69,10 +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),
-  fCutDCAToVertex(0),
-  fCutDCAToVertexXY(0),
+  fCutMaxDCAToVertexXY(0),
+  fCutMaxDCAToVertexZ(0),
+  fCutMinDCAToVertexXY(0),
+  fCutMinDCAToVertexZ(0),
+  fCutDCAToVertex2D(0),
   fPMin(0),
   fPMax(0),
   fPtMin(0),
@@ -107,11 +122,14 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   SetMaxCovDiagonalElements();                                     
   SetRequireTPCRefit();
   SetRequireITSRefit();
-  SetAcceptKingDaughters();
-  SetMinNsigmaToVertex();
-  SetRequireSigmaToVertex();
-  SetDCAToVertex();
-  SetDCAToVertexXY();
+  SetRequireITSStandAlone(kFALSE);
+  SetAcceptKinkDaughters();
+  SetMaxNsigmaToVertex();
+  SetMaxDCAToVertexXY();
+  SetMaxDCAToVertexZ();
+  SetDCAToVertex2D();
+  SetMinDCAToVertexXY();
+  SetMinDCAToVertexZ();
   SetPRange();
   SetPtRange();
   SetPxRange();
@@ -119,6 +137,9 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   SetPzRange();
   SetEtaRange();
   SetRapRange();
+  SetClusterRequirementITS(kSPD);
+  SetClusterRequirementITS(kSDD);
+  SetClusterRequirementITS(kSSD);
 
   SetHistogramsOn();
 }
@@ -137,10 +158,14 @@ AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
   fCutAcceptKinkDaughters(0),
   fCutRequireTPCRefit(0),
   fCutRequireITSRefit(0),
+  fCutRequireITSStandAlone(0),
   fCutNsigmaToVertex(0),
   fCutSigmaToVertexRequired(0),
-  fCutDCAToVertex(0),
-  fCutDCAToVertexXY(0),
+  fCutMaxDCAToVertexXY(0),
+  fCutMaxDCAToVertexZ(0),
+  fCutMinDCAToVertexXY(0),
+  fCutMinDCAToVertexZ(0),
+  fCutDCAToVertex2D(0),
   fPMin(0),
   fPMax(0),
   fPtMin(0),
@@ -237,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;
@@ -247,12 +275,17 @@ void AliESDtrackCuts::Init()
   fCutAcceptKinkDaughters = 0;
   fCutRequireTPCRefit = 0;
   fCutRequireITSRefit = 0;
+  fCutRequireITSStandAlone = 0;
 
   fCutNsigmaToVertex = 0;
   fCutSigmaToVertexRequired = 0;
-  fCutDCAToVertex = 0;
-  fCutDCAToVertexXY = 0;
+  fCutMaxDCAToVertexXY = 0;
+  fCutMaxDCAToVertexZ = 0;
+  fCutDCAToVertex2D = 0;
+  fCutMinDCAToVertexXY = 0;
+  fCutMinDCAToVertexZ = 0;
 
+  
   fPMin = 0;
   fPMax = 0;
   fPtMin = 0;
@@ -331,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;
@@ -340,11 +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.fCutDCAToVertex = fCutDCAToVertex;
-  target.fCutDCAToVertexXY = fCutDCAToVertexXY;
+  target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
+  target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
+  target.fCutDCAToVertex2D = fCutDCAToVertex2D;
+  target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
+  target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
 
   target.fPMin = fPMin;
   target.fPMax = fPMax;
@@ -470,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]);
@@ -484,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)
@@ -492,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)
@@ -522,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
   //
@@ -548,7 +589,6 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   // fTracks.fP   //GetMass
   // fTracks.fKinkIndexes
 
-
   UInt_t status = esdTrack->GetStatus();
 
   // getting quality parameters from the ESD track
@@ -574,10 +614,19 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
     AliDebug(1, "Estimated b resolution lower or equal zero!");
     bCov[0]=0; bCov[2]=0;
   }
-  Float_t dcaToVertex   = TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
+
   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];
@@ -644,18 +693,32 @@ 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 (dcaToVertex > fCutDCAToVertex)
+  if (fCutDCAToVertex2D && dcaToVertex > 1)
     cuts[21] = kTRUE;
-  if (dcaToVertexXY > fCutDCAToVertexXY)
+  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;}
 
 
 
@@ -729,6 +792,26 @@ AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
   return kTRUE;
 }
 
+//____________________________________________________________________
+Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
+{
+  // 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 kFALSE;
+}
+
 //____________________________________________________________________
 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
 {
@@ -744,6 +827,9 @@ AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
   if (!esd->GetPrimaryVertexTPC())
     return 0; // No TPC vertex no TPC tracks
 
+  if(!esd->GetPrimaryVertexTPC()->GetStatus())
+    return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
+
   AliESDtrack* track = esd->GetTrack(iTrack);
   if (!track)
     return 0;
@@ -788,6 +874,7 @@ TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); 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 = GetTPCOnlyTrack(esd, iTrack);
       if (!tpcTrack)
@@ -854,7 +941,7 @@ 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);
@@ -869,8 +956,8 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
 
     fhC11[i]                 = new TH1F("covMatrixDiagonal11","",2000,0,20);
     fhC22[i]                 = new TH1F("covMatrixDiagonal22","",2000,0,20);
-    fhC33[i]                 = new TH1F("covMatrixDiagonal33","",1000,0,1);
-    fhC44[i]                 = new TH1F("covMatrixDiagonal44","",1000,0,5);
+    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);
@@ -884,7 +971,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
 
     fhNSigmaToVertex[i]      = new TH1F("nSigmaToVertex","",500,0,10);
 
-    fhPt[i]                  = new TH1F("pt"     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
+    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");
@@ -898,16 +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");
-    fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2) in cm");
-    fhDXYvsDZ[i]->SetXTitle("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);
@@ -1188,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();
+}