[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) {
    hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
    
    if (cut)
      hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
    
    for (Int_t i=0; i<fNCuts; i++) {
      if (cuts[i])
 	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
      
      for (Int_t j=i; j<fNCuts; 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);
 	}
      }
    }
    

    hNClustersITS[0]->Fill(nClustersITS);        
    hNClustersTPC[0]->Fill(nClustersTPC);        
    hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
    hChi2PerClusterTPC[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]);                                  
    
    hDZ[0]->Fill(b[1]);     
    hDXY[0]->Fill(b[0]);    
    hDXYvsDZ[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]);
    }
  }

  //########################################################################  
  // cut the track!
  if (cut) return kFALSE;

  //########################################################################  
  // filling histograms after cut
  if (fHistogramsOn) {
    hNClustersITS[1]->Fill(nClustersITS);        
    hNClustersTPC[1]->Fill(nClustersTPC);        
    hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
    hChi2PerClusterTPC[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]);                                  
    
    hDZ[1]->Fill(b[1]);     
    hDXY[1]->Fill(b[0]);    
    hDXYvsDZ[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]);
  }
  
  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

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

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

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

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

  hDXYNormalized       = new TH1F*[2];
  hDZNormalized        = new TH1F*[2];
  hDXYvsDZNormalized   = new TH2F*[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);

  }
}

//____________________________________________________________________
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");
 
  hCutStatistics->Write();
  hCutCorrelation->Write();

  for (Int_t i=0; i<2; i++) {
    if (i==0)
      gDirectory->cd("before_cuts");
    else
      gDirectory->cd("after_cuts");
    
    hNClustersITS[i]        ->Write();
    hNClustersTPC[i]        ->Write();
    hChi2PerClusterITS[i]   ->Write();
    hChi2PerClusterTPC[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();
    
    hDXYNormalized[i]       ->Write();
    hDZNormalized[i]        ->Write();
    hDXYvsDZNormalized[i]   ->Write();

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

  gDirectory->cd("../");
}
Commit	Line	Data
0b75bef2	1	#include "AliESDtrackCuts.h"
	2
	3	#include <Riostream.h>
	4
	5	//____________________________________________________________________
	6	ClassImp(AliESDtrackCuts);
	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) {
179	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
180
181	if (cut)
182	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
183
184	for (Int_t i=0; i<fNCuts; i++) {
185	if (cuts[i])
186	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
187
188	for (Int_t j=i; j<fNCuts; j++) {
189	if (cuts[i] && cuts[j]) {
190	Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
191	Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
192	hCutCorrelation->Fill(x,y);
193	}
194	}
195	}
196
197
198	hNClustersITS[0]->Fill(nClustersITS);
199	hNClustersTPC[0]->Fill(nClustersTPC);
200	hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
201	hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
202
203	hC11[0]->Fill(extCov[0]);
204	hC22[0]->Fill(extCov[2]);
205	hC33[0]->Fill(extCov[5]);
206	hC44[0]->Fill(extCov[9]);
207	hC55[0]->Fill(extCov[14]);
208
209	hDZ[0]->Fill(b[1]);
210	hDXY[0]->Fill(b[0]);
211	hDXYvsDZ[0]->Fill(b[1],b[0]);
212
213	if (bRes[0]!=0 && bRes[1]!=0) {
214	hDZNormalized[0]->Fill(b[1]/bRes[1]);
215	hDXYNormalized[0]->Fill(b[0]/bRes[0]);
216	hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
217	}
218	}
219
220	//########################################################################
221	// cut the track!
222	if (cut) return kFALSE;
223
224	//########################################################################
225	// filling histograms after cut
226	if (fHistogramsOn) {
227	hNClustersITS[1]->Fill(nClustersITS);
228	hNClustersTPC[1]->Fill(nClustersTPC);
229	hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
230	hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
231
232	hC11[1]->Fill(extCov[0]);
233	hC22[1]->Fill(extCov[2]);
234	hC33[1]->Fill(extCov[5]);
235	hC44[1]->Fill(extCov[9]);
236	hC55[1]->Fill(extCov[14]);
237
238	hDZ[1]->Fill(b[1]);
239	hDXY[1]->Fill(b[0]);
240	hDXYvsDZ[1]->Fill(b[1],b[0]);
241
242	hDZNormalized[1]->Fill(b[1]/bRes[1]);
243	hDXYNormalized[1]->Fill(b[0]/bRes[0]);
244	hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
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	//____________________________________________________________________
268	void
269	AliESDtrackCuts::DefineHistograms(Int_t color) {
270
271	fHistogramsOn=kTRUE;
272
273	//###################################################################################
274	// defining histograms
275
276	hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
277
278	hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
279	hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
280
281	hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
282
283	for (Int_t i=0; i<fNCuts; i++) {
284	hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
285	hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
286	hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
287	}
288
289	hCutStatistics ->SetLineColor(color);
290	hCutCorrelation ->SetLineColor(color);
291	hCutStatistics ->SetLineWidth(2);
292	hCutCorrelation ->SetLineWidth(2);
293
294
295	hNClustersITS = new TH1F*[2];
296	hNClustersTPC = new TH1F*[2];
297	hChi2PerClusterITS = new TH1F*[2];
298	hChi2PerClusterTPC = new TH1F*[2];
299
300	hC11 = new TH1F*[2];
301	hC22 = new TH1F*[2];
302	hC33 = new TH1F*[2];
303	hC44 = new TH1F*[2];
304	hC55 = new TH1F*[2];
305
306	hDXY = new TH1F*[2];
307	hDZ = new TH1F*[2];
308	hDXYvsDZ = new TH2F*[2];
309
310	hDXYNormalized = new TH1F*[2];
311	hDZNormalized = new TH1F*[2];
312	hDXYvsDZNormalized = new TH2F*[2];
313
314
315	Char_t str[256];
316	for (Int_t i=0; i<2; i++) {
317	if (i==0) sprintf(str," ");
318	else sprintf(str,"_cut");
319
320	hNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
321	hNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
322	hChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
323	hChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
324
325	hC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
326	hC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
327	hC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
328	hC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
329	hC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
330
331	hDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10);
332	hDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10);
333	hDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
334
335	hDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
336	hDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
337	hDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
338
339
340	hNClustersITS[i] ->SetXTitle("n ITS clusters");
341	hNClustersTPC[i] ->SetXTitle("n TPC clusters");
342	hChi2PerClusterITS[i] ->SetXTitle("#Chi^{2} per ITS cluster");
343	hChi2PerClusterTPC[i] ->SetXTitle("#Chi^{2} per TPC cluster");
344
345	hC11[i] ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
346	hC22[i] ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
347	hC33[i] ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
348	hC44[i] ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
349	hC55[i] ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
350
351	hDXY[i] ->SetXTitle("transverse impact parameter");
352	hDZ[i] ->SetXTitle("longitudinal impact parameter");
353	hDXYvsDZ[i] ->SetXTitle("longitudinal impact parameter");
354	hDXYvsDZ[i] ->SetYTitle("transverse impact parameter");
355
356	hDXYNormalized[i] ->SetXTitle("normalized trans impact par");
357	hDZNormalized[i] ->SetXTitle("normalized long impact par");
358	hDXYvsDZNormalized[i] ->SetXTitle("normalized long impact par");
359	hDXYvsDZNormalized[i] ->SetYTitle("normalized trans impact par");
360
361	hNClustersITS[i] ->SetLineColor(color); hNClustersITS[i] ->SetLineWidth(2);
362	hNClustersTPC[i] ->SetLineColor(color); hNClustersTPC[i] ->SetLineWidth(2);
363	hChi2PerClusterITS[i] ->SetLineColor(color); hChi2PerClusterITS[i] ->SetLineWidth(2);
364	hChi2PerClusterTPC[i] ->SetLineColor(color); hChi2PerClusterTPC[i] ->SetLineWidth(2);
365
366	hC11[i] ->SetLineColor(color); hC11[i] ->SetLineWidth(2);
367	hC22[i] ->SetLineColor(color); hC22[i] ->SetLineWidth(2);
368	hC33[i] ->SetLineColor(color); hC33[i] ->SetLineWidth(2);
369	hC44[i] ->SetLineColor(color); hC44[i] ->SetLineWidth(2);
370	hC55[i] ->SetLineColor(color); hC55[i] ->SetLineWidth(2);
371
372	hDXY[i] ->SetLineColor(color); hDXY[i] ->SetLineWidth(2);
373	hDZ[i] ->SetLineColor(color); hDZ[i] ->SetLineWidth(2);
374
375	hDXYNormalized[i] ->SetLineColor(color); hDXYNormalized[i] ->SetLineWidth(2);
376	hDZNormalized[i] ->SetLineColor(color); hDZNormalized[i] ->SetLineWidth(2);
377
378	}
379	}
380
381	//____________________________________________________________________
382	void
49dc84d9	383	AliESDtrackCuts::Print(const Option_t*) const {
0b75bef2	384
	385	AliInfo("AliESDtrackCuts...");
	386	}
	387
	388
	389	//____________________________________________________________________
	390	void
	391	AliESDtrackCuts::SaveHistograms(Char_t* dir) {
	392
	393	if (!fHistogramsOn) {
	394	AliDebug(0, "Histograms not on - cannot save histograms!!!");
	395	return;
	396	}
	397
	398	gDirectory->mkdir(dir);
	399	gDirectory->cd(dir);
	400
	401	gDirectory->mkdir("before_cuts");
	402	gDirectory->mkdir("after_cuts");
	403
	404	hCutStatistics->Write();
	405	hCutCorrelation->Write();
	406
	407	for (Int_t i=0; i<2; i++) {
	408	if (i==0)
	409	gDirectory->cd("before_cuts");
	410	else
	411	gDirectory->cd("after_cuts");
	412
	413	hNClustersITS[i] ->Write();
	414	hNClustersTPC[i] ->Write();
	415	hChi2PerClusterITS[i] ->Write();
	416	hChi2PerClusterTPC[i] ->Write();
	417
	418	hC11[i] ->Write();
	419	hC22[i] ->Write();
	420	hC33[i] ->Write();
	421	hC44[i] ->Write();
	422	hC55[i] ->Write();
	423
	424	hDXY[i] ->Write();
	425	hDZ[i] ->Write();
	426	hDXYvsDZ[i] ->Write();
	427
	428	hDXYNormalized[i] ->Write();
	429	hDZNormalized[i] ->Write();
	430	hDXYvsDZNormalized[i] ->Write();
	431
	432	gDirectory->cd("../");
	433	}
	434
	435	gDirectory->cd("../");
	436	}
	437
	438
	439