[u/mrichter/AliRoot.git] / PWG0 / esdTrackCuts / AliESDtrackCuts.cxx

#include "AliESDtrackCuts.h"

#include <Riostream.h>

//____________________________________________________________________
ClassImp(AliESDtrackCuts)

//____________________________________________________________________
AliESDtrackCuts::AliESDtrackCuts() {

  //##############################################################################
  // setting default cuts

  SetMinNClustersTPC();
  SetMinNClustersITS();	    
  SetMaxChi2PerClusterTPC();
  SetMaxChi2PerClusterITS();  				    
  SetMaxCovDiagonalElements();  				    
  SetRequireTPCRefit();
  SetRequireITSRefit();
  SetAcceptKingDaughters();
  SetMinNsigmaToVertex();
  SetRequireSigmaToVertex();
  SetPRange();
  SetPtRange();
  SetPxRange();
  SetPyRange();
  SetPzRange();
  SetEtaRange();
  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";

}

//____________________________________________________________________
Bool_t 
AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
  // 
  // figure out if the tracks survives all the track cuts defined
  //

  UInt_t status = esdTrack->GetStatus();
  
  // 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)
    chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
  if (nClustersTPC!=0)
    chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);  

  Double_t extCov[15];
  esdTrack->GetExternalCovariance(extCov);  

  // getting the track to vertex parameters
  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 zero!");
    bCov[0]=0; bCov[1]=0;
  }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);

  // FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //
  // this is not correct - it will not give n sigma!!!
  // 
  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));  

  // getting the kinematic variables of the track 
  // (assuming the mass is known)
  Double_t p[3];
  esdTrack->GetPxPyPz(p);
  Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
  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.;
  if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
    eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
  if((energy != TMath::Abs(p[2]))&&(momentum != 0))
    y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));

  
  //########################################################################
  // cut the track?
  
  Bool_t cuts[fNCuts];
  for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
  
  // track quality cuts
  if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
    cuts[0]=kTRUE;
  if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
    cuts[1]=kTRUE;
  if (nClustersTPC<fCut_MinNClusterTPC) 
    cuts[2]=kTRUE;
  if (nClustersITS<fCut_MinNClusterITS) 
    cuts[3]=kTRUE;
  if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC) 
    cuts[4]=kTRUE; 
  if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS) 
    cuts[5]=kTRUE;
  if (extCov[0]  > fCut_MaxC11) 
    cuts[6]=kTRUE;  
  if (extCov[2]  > fCut_MaxC22) 
    cuts[7]=kTRUE;  
  if (extCov[5]  > fCut_MaxC33) 
    cuts[8]=kTRUE;  
  if (extCov[9]  > fCut_MaxC44) 
    cuts[9]=kTRUE;  
  if (extCov[14]  > fCut_MaxC55) 
    cuts[10]=kTRUE;  
  if (nSigmaToVertex > fCut_NsigmaToVertex) 
    cuts[11] = kTRUE;
  // if n sigma could not be calculated
  if (nSigmaToVertex<0 && fCut_SigmaToVertexRequired)   
    cuts[12]=kTRUE;
  if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
    cuts[13]=kTRUE;
  // track kinematics cut
  if((momentum < fPMin) || (momentum > fPMax)) 
    cuts[14]=kTRUE;
  if((pt < fPtMin) || (pt > fPtMax)) 
    cuts[15] = kTRUE;
  if((p[0] < fPxMin) || (p[0] > fPxMax)) 
    cuts[16] = kTRUE;
  if((p[1] < fPyMin) || (p[1] > fPyMax)) 
    cuts[17] = kTRUE;
  if((p[2] < fPzMin) || (p[2] > fPzMax)) 
    cuts[18] = kTRUE;
  if((eta < fEtaMin) || (eta > fEtaMax)) 
    cuts[19] = kTRUE;
  if((y < fRapMin) || (y > fRapMax)) 
    cuts[20] = kTRUE;

  Bool_t cut=kFALSE;
  for (Int_t i=0; i<fNCuts; i++) 
    if (cuts[i]) cut = kTRUE;
  
  //########################################################################
  // filling histograms
  if (fHistogramsOn) {
    fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
    
    if (cut)
      fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
    
    for (Int_t i=0; i<fNCuts; i++) {
      if (cuts[i])
 	fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
      
      for (Int_t j=i; j<fNCuts; j++) {
 	if (cuts[i] && cuts[j]) {
 	  Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
 	  Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
 	  fhCutCorrelation->Fill(x,y);
 	}
      }
    }
    

    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]);                                  
    
    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]);
    }
  }

  //########################################################################  
  // 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]);                                  
    
    fhDZ[1]->Fill(b[1]);     
    fhDXY[1]->Fill(b[0]);    
    fhDXYvsDZ[1]->Fill(b[1],b[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;
}

//____________________________________________________________________
TObjArray*
AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) {
  
  // returns an array of all tracks that pass the cuts
  fAcceptedTracks->Clear();
  
  // loop over esd tracks
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
    AliESDtrack* track = esd->GetTrack(iTrack);
    
    if(AcceptTrack(track)) fAcceptedTracks->Add(track);
  }

  return fAcceptedTracks;
}

//____________________________________________________________________
 void AliESDtrackCuts::DefineHistograms(Int_t color) {

   fHistogramsOn=kTRUE;

//   //###################################################################################
   // defining histograms

   fhCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);

   fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
   fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");

   fhCutCorrelation = 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++) {
     fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
     fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
     fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
   } 

  fhCutStatistics  ->SetLineColor(color);
  fhCutCorrelation ->SetLineColor(color);
  fhCutStatistics  ->SetLineWidth(2);
  fhCutCorrelation ->SetLineWidth(2);


  TH1F *fhNClustersITS  = new TH1F[2];
  TH1F *fhNClustersTPC  = new TH1F[2];
  TH1F *fhChi2PerClusterITS   = new TH1F[2];
  TH1F *fhChi2PerClusterTPC   = new TH1F[2];
  		       
  TH1F *fhC11                 = new TH1F[2];
  TH1F *fhC22                 = new TH1F[2];
  TH1F *fhC33                 = new TH1F[2];
  TH1F *fhC44                 = new TH1F[2];
  TH1F *fhC55                 = new TH1F[2];
  
  TH1F *fhDXY                 = new TH1F[2];
  TH1F *fhDZ		       = new TH1F[2];
  TH2F *fhDXYvsDZ	       = new TH2F[2];

  TH1F *fhDXYNormalized       = new TH1F[2];
  TH1F *fhDZNormalized        = new TH1F[2];
  TH2F *fhDXYvsDZNormalized   = new TH2F[2];


  Char_t str[256];
  for (Int_t i=0; i<2; i++) {
    if (i==0) sprintf(str," ");
    else sprintf(str,"_cut");

    fhNClustersITS[i]        = TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
    fhNClustersTPC[i]        = TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
    fhChi2PerClusterITS[i]   = TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
    fhChi2PerClusterTPC[i]   = TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
    
    fhC11[i]                 = TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
    fhC22[i]                 =  TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
    fhC33[i]                 =  TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
    fhC44[i]                 =  TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
    fhC55[i]                 =  TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
    
    fhDXY[i]                 =  TH1F(Form("dXY%s",str),"",500,-10,10);
    fhDZ[i]                  =  TH1F(Form("dZ%s",str),"",500,-10,10);
    fhDXYvsDZ[i]             =  TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);

    fhDXYNormalized[i]       =  TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
    fhDZNormalized[i]        =  TH1F(Form("dZNormalized%s",str),"",500,-10,10);
    fhDXYvsDZNormalized[i]   =  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(const Option_t*) const {

  AliInfo("AliESDtrackCuts...");
}


//____________________________________________________________________
void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
  
  if (!fHistogramsOn) {
    AliDebug(0, "Histograms not on - cannot save histograms!!!");
    return;
  }

  gDirectory->mkdir(dir);
  gDirectory->cd(dir);

  gDirectory->mkdir("before_cuts");
  gDirectory->mkdir("after_cuts");
 
  fhCutStatistics->Write();
  fhCutCorrelation->Write();

  for (Int_t i=0; i<2; i++) {
    if (i==0)
      gDirectory->cd("before_cuts");
    else
      gDirectory->cd("after_cuts");
    
    fhNClustersITS[i]        ->Write();
    fhNClustersTPC[i]        ->Write();
    fhChi2PerClusterITS[i]   ->Write();
    fhChi2PerClusterTPC[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();
    
    fhDXYNormalized[i]       ->Write();
    fhDZNormalized[i]        ->Write();
    fhDXYvsDZNormalized[i]   ->Write();

    gDirectory->cd("../");
  }

  gDirectory->cd("../");
}
Commit	Line	Data
0b75bef2	1	#include "AliESDtrackCuts.h"
	2
	3	#include <Riostream.h>
	4
	5	//____________________________________________________________________
f58f1a93	6	ClassImp(AliESDtrackCuts)
0b75bef2	7
	8	//____________________________________________________________________
	9	AliESDtrackCuts::AliESDtrackCuts() {
	10
	11	//##############################################################################
	12	// setting default cuts
	13
	14	SetMinNClustersTPC();
	15	SetMinNClustersITS();
	16	SetMaxChi2PerClusterTPC();
	17	SetMaxChi2PerClusterITS();
	18	SetMaxCovDiagonalElements();
	19	SetRequireTPCRefit();
	20	SetRequireITSRefit();
	21	SetAcceptKingDaughters();
	22	SetMinNsigmaToVertex();
	23	SetRequireSigmaToVertex();
	24	SetPRange();
	25	SetPtRange();
	26	SetPxRange();
	27	SetPyRange();
	28	SetPzRange();
	29	SetEtaRange();
	30	SetRapRange();
	31
	32	SetHistogramsOn();
	33
	34	// set the cut names
	35	fCutNames[0] = "require TPC refit";
	36	fCutNames[1] = "require ITS refit";
	37	fCutNames[2] = "n clusters TPC";
	38	fCutNames[3] = "n clusters ITS";
	39	fCutNames[4] = "#Chi^{2}/clusters TPC";
	40	fCutNames[5] = "#Chi^{2}/clusters ITS";
	41	fCutNames[6] = "cov 11";
	42	fCutNames[7] = "cov 22";
	43	fCutNames[8] = "cov 33";
	44	fCutNames[9] = "cov 44";
	45	fCutNames[10] = "cov 55";
	46	fCutNames[11] = "trk-to-vtx";
	47	fCutNames[12] = "trk-to-vtx failed";
	48	fCutNames[13] = "kink daughters";
	49
	50	fCutNames[14] = "p";
	51	fCutNames[15] = "p_{T}";
	52	fCutNames[16] = "p_{x}";
	53	fCutNames[17] = "p_{y}";
	54	fCutNames[18] = "p_{z}";
	55	fCutNames[19] = "y";
	56	fCutNames[20] = "eta";
	57
	58	}
	59
0b75bef2	60	//____________________________________________________________________
	61	Bool_t
	62	AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) {
	63	//
	64	// figure out if the tracks survives all the track cuts defined
	65	//
	66
	67	UInt_t status = esdTrack->GetStatus();
	68
	69	// getting quality parameters from the ESD track
	70	Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
	71	Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
	72
	73	Float_t chi2PerClusterITS = -1;
	74	Float_t chi2PerClusterTPC = -1;
	75	if (nClustersITS!=0)
	76	chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
	77	if (nClustersTPC!=0)
	78	chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
	79
	80	Double_t extCov[15];
	81	esdTrack->GetExternalCovariance(extCov);
	82
	83	// getting the track to vertex parameters
	84	Float_t b[2];
	85	Float_t bRes[2];
	86	Float_t bCov[3];
	87	esdTrack->GetImpactParameters(b,bCov);
	88	if (bCov[0]<=0 \|\| bCov[2]<=0) {
	89	AliDebug(1, "Estimated b resolution zero!");
	90	bCov[0]=0; bCov[1]=0;
	91	}
	92	bRes[0] = TMath::Sqrt(bCov[0]);
	93	bRes[1] = TMath::Sqrt(bCov[2]);
	94
	95	// FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	96	//
	97	// this is not correct - it will not give n sigma!!!
	98	//
	99	Float_t nSigmaToVertex = -1;
	100	if (bRes[0]!=0 && bRes[1]!=0)
	101	nSigmaToVertex = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	102
	103	// getting the kinematic variables of the track
	104	// (assuming the mass is known)
	105	Double_t p[3];
	106	esdTrack->GetPxPyPz(p);
	107	Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
	108	Float_t pt = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
	109	Float_t energy = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
	110
	111	//y-eta related calculations
	112	Float_t eta = -100.;
	113	Float_t y = -100.;
	114	if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
	115	eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
	116	if((energy != TMath::Abs(p[2]))&&(momentum != 0))
	117	y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
	118
	119
	120	//########################################################################
	121	// cut the track?
	122
	123	Bool_t cuts[fNCuts];
124	for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
125
126	// track quality cuts
127	if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
128	cuts[0]=kTRUE;
129	if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
130	cuts[1]=kTRUE;
131	if (nClustersTPC<fCut_MinNClusterTPC)
132	cuts[2]=kTRUE;
133	if (nClustersITS<fCut_MinNClusterITS)
134	cuts[3]=kTRUE;
135	if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
136	cuts[4]=kTRUE;
137	if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
138	cuts[5]=kTRUE;
139	if (extCov[0] > fCut_MaxC11)
140	cuts[6]=kTRUE;
141	if (extCov[2] > fCut_MaxC22)
142	cuts[7]=kTRUE;
143	if (extCov[5] > fCut_MaxC33)
144	cuts[8]=kTRUE;
145	if (extCov[9] > fCut_MaxC44)
146	cuts[9]=kTRUE;
147	if (extCov[14] > fCut_MaxC55)
148	cuts[10]=kTRUE;
149	if (nSigmaToVertex > fCut_NsigmaToVertex)
150	cuts[11] = kTRUE;
151	// if n sigma could not be calculated
152	if (nSigmaToVertex<0 && fCut_SigmaToVertexRequired)
153	cuts[12]=kTRUE;
154	if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
155	cuts[13]=kTRUE;
156	// track kinematics cut
157	if((momentum < fPMin) \|\| (momentum > fPMax))
158	cuts[14]=kTRUE;
159	if((pt < fPtMin) \|\| (pt > fPtMax))
160	cuts[15] = kTRUE;
161	if((p[0] < fPxMin) \|\| (p[0] > fPxMax))
162	cuts[16] = kTRUE;
163	if((p[1] < fPyMin) \|\| (p[1] > fPyMax))
164	cuts[17] = kTRUE;
165	if((p[2] < fPzMin) \|\| (p[2] > fPzMax))
166	cuts[18] = kTRUE;
167	if((eta < fEtaMin) \|\| (eta > fEtaMax))
168	cuts[19] = kTRUE;
169	if((y < fRapMin) \|\| (y > fRapMax))
170	cuts[20] = kTRUE;
171
172	Bool_t cut=kFALSE;
173	for (Int_t i=0; i<fNCuts; i++)
174	if (cuts[i]) cut = kTRUE;
175
176	//########################################################################
177	// filling histograms
178	if (fHistogramsOn) {
f58f1a93	179	fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
0b75bef2	180
0b75bef2	181	if (cut)
f58f1a93	182	fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
0b75bef2	183
	184	for (Int_t i=0; i<fNCuts; i++) {
	185	if (cuts[i])
f58f1a93	186	fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
0b75bef2	187
	188	for (Int_t j=i; j<fNCuts; j++) {
	189	if (cuts[i] && cuts[j]) {
f58f1a93	190	Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
	191	Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
	192	fhCutCorrelation->Fill(x,y);
0b75bef2	193	}
	194	}
	195	}
	196
	197
f58f1a93	198	fhNClustersITS[0]->Fill(nClustersITS);
	199	fhNClustersTPC[0]->Fill(nClustersTPC);
	200	fhChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
	201	fhChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
0b75bef2	202
f58f1a93	203	fhC11[0]->Fill(extCov[0]);
	204	fhC22[0]->Fill(extCov[2]);
	205	fhC33[0]->Fill(extCov[5]);
	206	fhC44[0]->Fill(extCov[9]);
	207	fhC55[0]->Fill(extCov[14]);
0b75bef2	208
f58f1a93	209	fhDZ[0]->Fill(b[1]);
	210	fhDXY[0]->Fill(b[0]);
	211	fhDXYvsDZ[0]->Fill(b[1],b[0]);
0b75bef2	212
0b75bef2	213	if (bRes[0]!=0 && bRes[1]!=0) {
f58f1a93	214	fhDZNormalized[0]->Fill(b[1]/bRes[1]);
	215	fhDXYNormalized[0]->Fill(b[0]/bRes[0]);
	216	fhDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
0b75bef2	217	}
	218	}
	219
	220	//########################################################################
	221	// cut the track!
	222	if (cut) return kFALSE;
	223
	224	//########################################################################
	225	// filling histograms after cut
	226	if (fHistogramsOn) {
f58f1a93	227	fhNClustersITS[1]->Fill(nClustersITS);
	228	fhNClustersTPC[1]->Fill(nClustersTPC);
	229	fhChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
	230	fhChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
0b75bef2	231
f58f1a93	232	fhC11[1]->Fill(extCov[0]);
	233	fhC22[1]->Fill(extCov[2]);
	234	fhC33[1]->Fill(extCov[5]);
	235	fhC44[1]->Fill(extCov[9]);
	236	fhC55[1]->Fill(extCov[14]);
0b75bef2	237
f58f1a93	238	fhDZ[1]->Fill(b[1]);
	239	fhDXY[1]->Fill(b[0]);
	240	fhDXYvsDZ[1]->Fill(b[1],b[0]);
0b75bef2	241
f58f1a93	242	fhDZNormalized[1]->Fill(b[1]/bRes[1]);
	243	fhDXYNormalized[1]->Fill(b[0]/bRes[0]);
	244	fhDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
0b75bef2	245	}
	246
	247	return kTRUE;
	248	}
	249
	250	//____________________________________________________________________
	251	TObjArray*
	252	AliESDtrackCuts::GetAcceptedTracks(AliESD* esd) {
	253
	254	// returns an array of all tracks that pass the cuts
	255	fAcceptedTracks->Clear();
	256
	257	// loop over esd tracks
	258	for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
	259	AliESDtrack* track = esd->GetTrack(iTrack);
	260
	261	if(AcceptTrack(track)) fAcceptedTracks->Add(track);
	262	}
	263
	264	return fAcceptedTracks;
	265	}
	266
	267	//____________________________________________________________________
f58f1a93	268	void AliESDtrackCuts::DefineHistograms(Int_t color) {
0b75bef2	269
f58f1a93	270	fHistogramsOn=kTRUE;
0b75bef2	271
f58f1a93	272	// //###################################################################################
f58f1a93	273	// defining histograms
0b75bef2	274
f58f1a93	275	fhCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
0b75bef2	276
f58f1a93	277	fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
f58f1a93	278	fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
0b75bef2	279
f58f1a93	280	fhCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
0b75bef2	281
f58f1a93	282	for (Int_t i=0; i<fNCuts; i++) {
	283	fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
	284	fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
	285	fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
	286	}
	287
	288	fhCutStatistics ->SetLineColor(color);
	289	fhCutCorrelation ->SetLineColor(color);
	290	fhCutStatistics ->SetLineWidth(2);
	291	fhCutCorrelation ->SetLineWidth(2);
	292
	293
	294	TH1F *fhNClustersITS = new TH1F[2];
	295	TH1F *fhNClustersTPC = new TH1F[2];
	296	TH1F *fhChi2PerClusterITS = new TH1F[2];
	297	TH1F *fhChi2PerClusterTPC = new TH1F[2];
0b75bef2	298
f58f1a93	299	TH1F *fhC11 = new TH1F[2];
	300	TH1F *fhC22 = new TH1F[2];
	301	TH1F *fhC33 = new TH1F[2];
	302	TH1F *fhC44 = new TH1F[2];
	303	TH1F *fhC55 = new TH1F[2];
0b75bef2	304
f58f1a93	305	TH1F *fhDXY = new TH1F[2];
	306	TH1F *fhDZ = new TH1F[2];
	307	TH2F *fhDXYvsDZ = new TH2F[2];
0b75bef2	308
f58f1a93	309	TH1F *fhDXYNormalized = new TH1F[2];
	310	TH1F *fhDZNormalized = new TH1F[2];
	311	TH2F *fhDXYvsDZNormalized = new TH2F[2];
0b75bef2	312
	313
	314	Char_t str[256];
	315	for (Int_t i=0; i<2; i++) {
	316	if (i==0) sprintf(str," ");
	317	else sprintf(str,"_cut");
	318
f58f1a93	319	fhNClustersITS[i] = TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
	320	fhNClustersTPC[i] = TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
	321	fhChi2PerClusterITS[i] = TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
	322	fhChi2PerClusterTPC[i] = TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
0b75bef2	323
f58f1a93	324	fhC11[i] = TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
	325	fhC22[i] = TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
	326	fhC33[i] = TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
	327	fhC44[i] = TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
	328	fhC55[i] = TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
0b75bef2	329
f58f1a93	330	fhDXY[i] = TH1F(Form("dXY%s",str),"",500,-10,10);
	331	fhDZ[i] = TH1F(Form("dZ%s",str),"",500,-10,10);
	332	fhDXYvsDZ[i] = TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
0b75bef2	333
f58f1a93	334	fhDXYNormalized[i] = TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
	335	fhDZNormalized[i] = TH1F(Form("dZNormalized%s",str),"",500,-10,10);
	336	fhDXYvsDZNormalized[i] = TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
0b75bef2	337
0b75bef2	338
f58f1a93	339	fhNClustersITS[i].SetXTitle("n ITS clusters");
	340	fhNClustersTPC[i].SetXTitle("n TPC clusters");
	341	fhChi2PerClusterITS[i].SetXTitle("#Chi^{2} per ITS cluster");
	342	fhChi2PerClusterTPC[i].SetXTitle("#Chi^{2} per TPC cluster");
0b75bef2	343
f58f1a93	344	fhC11[i].SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
	345	fhC22[i].SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
	346	fhC33[i].SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
	347	fhC44[i].SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
	348	fhC55[i].SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
0b75bef2	349
f58f1a93	350	fhDXY[i].SetXTitle("transverse impact parameter");
	351	fhDZ[i].SetXTitle("longitudinal impact parameter");
	352	fhDXYvsDZ[i].SetXTitle("longitudinal impact parameter");
	353	fhDXYvsDZ[i].SetYTitle("transverse impact parameter");
	354
	355	fhDXYNormalized[i].SetXTitle("normalized trans impact par");
	356	fhDZNormalized[i].SetXTitle("normalized long impact par");
	357	fhDXYvsDZNormalized[i].SetXTitle("normalized long impact par");
	358	fhDXYvsDZNormalized[i].SetYTitle("normalized trans impact par");
	359
	360	fhNClustersITS[i].SetLineColor(color); fhNClustersITS[i].SetLineWidth(2);
	361	fhNClustersTPC[i].SetLineColor(color); fhNClustersTPC[i].SetLineWidth(2);
	362	fhChi2PerClusterITS[i].SetLineColor(color); fhChi2PerClusterITS[i].SetLineWidth(2);
	363	fhChi2PerClusterTPC[i].SetLineColor(color); fhChi2PerClusterTPC[i].SetLineWidth(2);
0b75bef2	364
f58f1a93	365	fhC11[i].SetLineColor(color); fhC11[i].SetLineWidth(2);
	366	fhC22[i].SetLineColor(color); fhC22[i].SetLineWidth(2);
	367	fhC33[i].SetLineColor(color); fhC33[i].SetLineWidth(2);
	368	fhC44[i].SetLineColor(color); fhC44[i].SetLineWidth(2);
	369	fhC55[i].SetLineColor(color); fhC55[i].SetLineWidth(2);
0b75bef2	370
f58f1a93	371	fhDXY[i].SetLineColor(color); fhDXY[i].SetLineWidth(2);
f58f1a93	372	fhDZ[i].SetLineColor(color); fhDZ[i].SetLineWidth(2);
0b75bef2	373
f58f1a93	374	fhDXYNormalized[i].SetLineColor(color); fhDXYNormalized[i].SetLineWidth(2);
f58f1a93	375	fhDZNormalized[i].SetLineColor(color); fhDZNormalized[i].SetLineWidth(2);
0b75bef2	376
	377	}
	378	}
	379
	380	//____________________________________________________________________
	381	void
49dc84d9	382	AliESDtrackCuts::Print(const Option_t*) const {
0b75bef2	383
	384	AliInfo("AliESDtrackCuts...");
	385	}
	386
	387
	388	//____________________________________________________________________
f58f1a93	389	void AliESDtrackCuts::SaveHistograms(Char_t* dir) {
0b75bef2	390
	391	if (!fHistogramsOn) {
	392	AliDebug(0, "Histograms not on - cannot save histograms!!!");
	393	return;
	394	}
	395
	396	gDirectory->mkdir(dir);
	397	gDirectory->cd(dir);
	398
	399	gDirectory->mkdir("before_cuts");
	400	gDirectory->mkdir("after_cuts");
	401
f58f1a93	402	fhCutStatistics->Write();
f58f1a93	403	fhCutCorrelation->Write();
0b75bef2	404
	405	for (Int_t i=0; i<2; i++) {
	406	if (i==0)
	407	gDirectory->cd("before_cuts");
	408	else
	409	gDirectory->cd("after_cuts");
	410
f58f1a93	411	fhNClustersITS[i] ->Write();
	412	fhNClustersTPC[i] ->Write();
	413	fhChi2PerClusterITS[i] ->Write();
	414	fhChi2PerClusterTPC[i] ->Write();
0b75bef2	415
f58f1a93	416	fhC11[i] ->Write();
	417	fhC22[i] ->Write();
	418	fhC33[i] ->Write();
	419	fhC44[i] ->Write();
	420	fhC55[i] ->Write();
	421
	422	fhDXY[i] ->Write();
	423	fhDZ[i] ->Write();
	424	fhDXYvsDZ[i] ->Write();
0b75bef2	425
f58f1a93	426	fhDXYNormalized[i] ->Write();
	427	fhDZNormalized[i] ->Write();
	428	fhDXYvsDZNormalized[i] ->Write();
0b75bef2	429
	430	gDirectory->cd("../");
	431	}
	432
	433	gDirectory->cd("../");
	434	}
	435
	436
	437