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

#include "ESDtrackQualityCuts.h"

#include <Riostream.h>

//____________________________________________________________________
ClassImp(ESDtrackQualityCuts);

//____________________________________________________________________
ESDtrackQualityCuts::ESDtrackQualityCuts() {

  //##############################################################################
  // setting default cuts
  SetMinNClustersTPC();
  SetMinNClustersITS();	    
  SetMaxChi2PerClusterTPC();
  SetMaxChi2PerClusterITS();  				    
  SetMaxCovDiagonalElements();  				    
  SetRequireTPCRefit();
  SetRequireITSRefit();
  SetAcceptKingDaughters();
  SetMinNsigmaToVertex();
  SetRequireSigmaToVertex();

  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";
}

//____________________________________________________________________
Bool_t 
ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {

  // relate the track to the new vertex
  esdTrack->RelateToVertex(esdVtx, field, 999999);
  
  return AcceptTrack(esdTrack);
}

//____________________________________________________________________
Bool_t 
ESDtrackQualityCuts::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);
} 

//____________________________________________________________________
Bool_t 
ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack) {

  UInt_t status = esdTrack->GetStatus();
  
  // getting variables from ESD
  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);  

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

  //########################################################################
  // cut the track?
  
  Bool_t cuts[fNCuts];
  for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;

  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==-1 && fCut_SigmaToVertexRequired)   
    cuts[12]=kTRUE;
  if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) 
    cuts[13]=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;
}

//____________________________________________________________________
void 
ESDtrackQualityCuts::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 
ESDtrackQualityCuts::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
75ec0f41	1	#include "ESDtrackQualityCuts.h"
	2
	3	#include <Riostream.h>
	4
	5	//____________________________________________________________________
	6	ClassImp(ESDtrackQualityCuts);
	7
	8	//____________________________________________________________________
	9	ESDtrackQualityCuts::ESDtrackQualityCuts() {
	10
	11	//##############################################################################
	12	// setting default cuts
	13	SetMinNClustersTPC();
	14	SetMinNClustersITS();
	15	SetMaxChi2PerClusterTPC();
	16	SetMaxChi2PerClusterITS();
	17	SetMaxCovDiagonalElements();
	18	SetRequireTPCRefit();
	19	SetRequireITSRefit();
	20	SetAcceptKingDaughters();
	21	SetMinNsigmaToVertex();
	22	SetRequireSigmaToVertex();
	23
	24	SetHistogramsOn();
	25
	26	// set the cut names
	27	fCutNames[0] = "require TPC refit";
	28	fCutNames[1] = "require ITS refit";
	29	fCutNames[2] = "n clusters TPC";
	30	fCutNames[3] = "n clusters ITS";
	31	fCutNames[4] = "#Chi^{2}/clusters TPC";
	32	fCutNames[5] = "#Chi^{2}/clusters ITS";
	33	fCutNames[6] = "cov 11";
	34	fCutNames[7] = "cov 22";
	35	fCutNames[8] = "cov 33";
	36	fCutNames[9] = "cov 44";
	37	fCutNames[10] = "cov 55";
	38	fCutNames[11] = "trk-to-vtx";
	39	fCutNames[12] = "trk-to-vtx failed";
	40	fCutNames[13] = "kink daughters";
	41	}
	42
	43	//____________________________________________________________________
	44	Bool_t
	45	ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, AliESDVertex* esdVtx, Double_t field) {
	46
	47	// relate the track to the new vertex
	48	esdTrack->RelateToVertex(esdVtx, field, 999999);
	49
	50	return AcceptTrack(esdTrack);
	51	}
	52
	53	//____________________________________________________________________
	54	Bool_t
	55	ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack, Double_t* vtx, Double_t* vtx_res, Double_t field) {
	56
	57	AliESDVertex* esdVtx = new AliESDVertex(vtx, vtx_res,"new vertex");
	58	esdTrack->RelateToVertex(esdVtx, field, 1e99);
	59	return AcceptTrack(esdTrack);
	60	}
	61
	62	//____________________________________________________________________
	63	Bool_t
	64	ESDtrackQualityCuts::AcceptTrack(AliESDtrack* esdTrack) {
65
66	UInt_t status = esdTrack->GetStatus();
67
68	// getting variables from ESD
69	Int_t nClustersITS = esdTrack->GetITSclusters(fIdxInt);
70	Int_t nClustersTPC = esdTrack->GetTPCclusters(fIdxInt);
71
72	Float_t chi2PerClusterITS = -1;
73	Float_t chi2PerClusterTPC = -1;
74	if (nClustersITS!=0)
75	chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
76	if (nClustersTPC!=0)
77	chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
78
79	Double_t extCov[15];
80	esdTrack->GetExternalCovariance(extCov);
81
82	Float_t b[2];
83	Float_t bRes[2];
84	Float_t bCov[3];
85	esdTrack->GetImpactParameters(b,bCov);
86	if (bCov[0]<=0 \|\| bCov[2]<=0) {
87	AliDebug(1, "Estimated b resolution zero!");
88	bCov[0]=0; bCov[1]=0;
89	}
90	bRes[0] = TMath::Sqrt(bCov[0]);
91	bRes[1] = TMath::Sqrt(bCov[2]);
92
93	//
94	// FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95	//
96	// this is not correct - it will not give n sigma!!!
97	//
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
104	//########################################################################
105	// cut the track?
106
107	Bool_t cuts[fNCuts];
108	for (Int_t i=0; i<fNCuts; i++) cuts[i]=kFALSE;
109
110	if (fCut_RequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
111	cuts[0]=kTRUE;
112	if (fCut_RequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
113	cuts[1]=kTRUE;
114	if (nClustersTPC<fCut_MinNClusterTPC)
115	cuts[2]=kTRUE;
116	if (nClustersITS<fCut_MinNClusterITS)
117	cuts[3]=kTRUE;
118	if (chi2PerClusterTPC>fCut_MaxChi2PerClusterTPC)
119	cuts[4]=kTRUE;
120	if (chi2PerClusterITS>fCut_MaxChi2PerClusterITS)
121	cuts[5]=kTRUE;
122	if (extCov[0] > fCut_MaxC11)
123	cuts[6]=kTRUE;
124	if (extCov[2] > fCut_MaxC22)
125	cuts[7]=kTRUE;
126	if (extCov[5] > fCut_MaxC33)
127	cuts[8]=kTRUE;
128	if (extCov[9] > fCut_MaxC44)
129	cuts[9]=kTRUE;
130	if (extCov[14] > fCut_MaxC55)
131	cuts[10]=kTRUE;
132	if (nSigmaToVertex > fCut_NsigmaToVertex)
133	cuts[11] = kTRUE;
134	// if n sigma could not be calculated
135	if (nSigmaToVertex==-1 && fCut_SigmaToVertexRequired)
136	cuts[12]=kTRUE;
137	if (!fCut_AcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
138	cuts[13]=kTRUE;
139
140	Bool_t cut=kFALSE;
141	for (Int_t i=0; i<fNCuts; i++)
142	if (cuts[i]) cut = kTRUE;
143
144	//########################################################################
145	// filling histograms
146	if (fHistogramsOn) {
147	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n tracks")));
148
149	if (cut)
150	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin("n cut tracks")));
151
152	for (Int_t i=0; i<fNCuts; i++) {
153	if (cuts[i])
154	hCutStatistics->Fill(hCutStatistics->GetBinCenter(hCutStatistics->GetXaxis()->FindBin(fCutNames[i])));
155
156	for (Int_t j=i; j<fNCuts; j++) {
157	if (cuts[i] && cuts[j]) {
158	Float_t x = hCutCorrelation->GetXaxis()->GetBinCenter(hCutCorrelation->GetXaxis()->FindBin(fCutNames[i]));
159	Float_t y = hCutCorrelation->GetYaxis()->GetBinCenter(hCutCorrelation->GetYaxis()->FindBin(fCutNames[j]));
160	hCutCorrelation->Fill(x,y);
161	}
162	}
163	}
164
165
166	hNClustersITS[0]->Fill(nClustersITS);
167	hNClustersTPC[0]->Fill(nClustersTPC);
168	hChi2PerClusterITS[0]->Fill(chi2PerClusterITS);
169	hChi2PerClusterTPC[0]->Fill(chi2PerClusterTPC);
170
171	hC11[0]->Fill(extCov[0]);
172	hC22[0]->Fill(extCov[2]);
173	hC33[0]->Fill(extCov[5]);
174	hC44[0]->Fill(extCov[9]);
175	hC55[0]->Fill(extCov[14]);
176
177	hDZ[0]->Fill(b[1]);
178	hDXY[0]->Fill(b[0]);
179	hDXYvsDZ[0]->Fill(b[1],b[0]);
180
181	if (bRes[0]!=0 && bRes[1]!=0) {
182	hDZNormalized[0]->Fill(b[1]/bRes[1]);
183	hDXYNormalized[0]->Fill(b[0]/bRes[0]);
184	hDXYvsDZNormalized[0]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
185	}
186	}
187
188	//########################################################################
189	// cut the track!
190	if (cut) return kFALSE;
191
192	//########################################################################
193	// filling histograms after cut
194	if (fHistogramsOn) {
195	hNClustersITS[1]->Fill(nClustersITS);
196	hNClustersTPC[1]->Fill(nClustersTPC);
197	hChi2PerClusterITS[1]->Fill(chi2PerClusterITS);
198	hChi2PerClusterTPC[1]->Fill(chi2PerClusterTPC);
199
200	hC11[1]->Fill(extCov[0]);
201	hC22[1]->Fill(extCov[2]);
202	hC33[1]->Fill(extCov[5]);
203	hC44[1]->Fill(extCov[9]);
204	hC55[1]->Fill(extCov[14]);
205
206	hDZ[1]->Fill(b[1]);
207	hDXY[1]->Fill(b[0]);
208	hDXYvsDZ[1]->Fill(b[1],b[0]);
209
210	hDZNormalized[1]->Fill(b[1]/bRes[1]);
211	hDXYNormalized[1]->Fill(b[0]/bRes[0]);
212	hDXYvsDZNormalized[1]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
213	}
214
215	return kTRUE;
216	}
217
218	//____________________________________________________________________
219	void
220	ESDtrackQualityCuts::DefineHistograms(Int_t color) {
221
222	fHistogramsOn=kTRUE;
223
224	//###################################################################################
225	// defining histograms
226
227	hCutStatistics = new TH1F("cut_statistics","cut statistics",fNCuts+4,-0.5,fNCuts+3.5);
228
229	hCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
230	hCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
231
232	hCutCorrelation = new TH2F("cut_correlation","cut correlation",fNCuts,-0.5,fNCuts-0.5,fNCuts,-0.5,fNCuts-0.5);;
233
234	for (Int_t i=0; i<fNCuts; i++) {
235	hCutStatistics->GetXaxis()->SetBinLabel(i+4,fCutNames[i]);
236	hCutCorrelation->GetXaxis()->SetBinLabel(i+1,fCutNames[i]);
237	hCutCorrelation->GetYaxis()->SetBinLabel(i+1,fCutNames[i]);
238	}
239
240	hCutStatistics ->SetLineColor(color);
241	hCutCorrelation ->SetLineColor(color);
242	hCutStatistics ->SetLineWidth(2);
243	hCutCorrelation ->SetLineWidth(2);
244
245
246	hNClustersITS = new TH1F*[2];
247	hNClustersTPC = new TH1F*[2];
248	hChi2PerClusterITS = new TH1F*[2];
249	hChi2PerClusterTPC = new TH1F*[2];
250
251	hC11 = new TH1F*[2];
252	hC22 = new TH1F*[2];
253	hC33 = new TH1F*[2];
254	hC44 = new TH1F*[2];
255	hC55 = new TH1F*[2];
256
257	hDXY = new TH1F*[2];
258	hDZ = new TH1F*[2];
259	hDXYvsDZ = new TH2F*[2];
260
261	hDXYNormalized = new TH1F*[2];
262	hDZNormalized = new TH1F*[2];
263	hDXYvsDZNormalized = new TH2F*[2];
264
265
266
267	Char_t str[256];
268	for (Int_t i=0; i<2; i++) {
269	if (i==0) sprintf(str,"");
270	else sprintf(str,"_cut");
271
272	hNClustersITS[i] = new TH1F(Form("nClustersITS%s",str),"",8,-0.5,7.5);
273	hNClustersTPC[i] = new TH1F(Form("nClustersTPC%s",str),"",165,-0.5,164.5);
274	hChi2PerClusterITS[i] = new TH1F(Form("chi2PerClusterITS%s",str),"",500,0,10);
275	hChi2PerClusterTPC[i] = new TH1F(Form("chi2PerClusterTPC%s",str),"",500,0,10);
276
277	hC11[i] = new TH1F(Form("covMatrixDiagonal11%s",str),"",1000,0,5);
278	hC22[i] = new TH1F(Form("covMatrixDiagonal22%s",str),"",1000,0,5);
279	hC33[i] = new TH1F(Form("covMatrixDiagonal33%s",str),"",1000,0,0.5);
280	hC44[i] = new TH1F(Form("covMatrixDiagonal44%s",str),"",1000,0,5);
281	hC55[i] = new TH1F(Form("covMatrixDiagonal55%s",str),"",1000,0,5);
282
283	hDXY[i] = new TH1F(Form("dXY%s",str),"",500,-10,10);
284	hDZ[i] = new TH1F(Form("dZ%s",str),"",500,-10,10);
285	hDXYvsDZ[i] = new TH2F(Form("dXYvsDZ%s",str),"",200,-10,10,200,-10,10);
286
287	hDXYNormalized[i] = new TH1F(Form("dXYNormalized%s",str),"",500,-10,10);
288	hDZNormalized[i] = new TH1F(Form("dZNormalized%s",str),"",500,-10,10);
289	hDXYvsDZNormalized[i] = new TH2F(Form("dXYvsDZNormalized%s",str),"",200,-10,10,200,-10,10);
290
291
292	hNClustersITS[i] ->SetXTitle("n ITS clusters");
293	hNClustersTPC[i] ->SetXTitle("n TPC clusters");
294	hChi2PerClusterITS[i] ->SetXTitle("#Chi^{2} per ITS cluster");
295	hChi2PerClusterTPC[i] ->SetXTitle("#Chi^{2} per TPC cluster");
296
297	hC11[i] ->SetXTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
298	hC22[i] ->SetXTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
299	hC33[i] ->SetXTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
300	hC44[i] ->SetXTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
301	hC55[i] ->SetXTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
302
303	hDXY[i] ->SetXTitle("transverse impact parameter");
304	hDZ[i] ->SetXTitle("longitudinal impact parameter");
305	hDXYvsDZ[i] ->SetXTitle("longitudinal impact parameter");
306	hDXYvsDZ[i] ->SetYTitle("transverse impact parameter");
307
308	hDXYNormalized[i] ->SetXTitle("normalized trans impact par");
309	hDZNormalized[i] ->SetXTitle("normalized long impact par");
310	hDXYvsDZNormalized[i] ->SetXTitle("normalized long impact par");
311	hDXYvsDZNormalized[i] ->SetYTitle("normalized trans impact par");
312
313	hNClustersITS[i] ->SetLineColor(color); hNClustersITS[i] ->SetLineWidth(2);
314	hNClustersTPC[i] ->SetLineColor(color); hNClustersTPC[i] ->SetLineWidth(2);
315	hChi2PerClusterITS[i] ->SetLineColor(color); hChi2PerClusterITS[i] ->SetLineWidth(2);
316	hChi2PerClusterTPC[i] ->SetLineColor(color); hChi2PerClusterTPC[i] ->SetLineWidth(2);
317
318	hC11[i] ->SetLineColor(color); hC11[i] ->SetLineWidth(2);
319	hC22[i] ->SetLineColor(color); hC22[i] ->SetLineWidth(2);
320	hC33[i] ->SetLineColor(color); hC33[i] ->SetLineWidth(2);
321	hC44[i] ->SetLineColor(color); hC44[i] ->SetLineWidth(2);
322	hC55[i] ->SetLineColor(color); hC55[i] ->SetLineWidth(2);
323
324	hDXY[i] ->SetLineColor(color); hDXY[i] ->SetLineWidth(2);
325	hDZ[i] ->SetLineColor(color); hDZ[i] ->SetLineWidth(2);
326
327	hDXYNormalized[i] ->SetLineColor(color); hDXYNormalized[i] ->SetLineWidth(2);
328	hDZNormalized[i] ->SetLineColor(color); hDZNormalized[i] ->SetLineWidth(2);
329
330	}
331	}
332
333	//____________________________________________________________________
334	void
335	ESDtrackQualityCuts::SaveHistograms(Char_t* dir) {
336
337	if (!fHistogramsOn) {
338	AliDebug(0, "Histograms not on - cannot save histograms!!!");
339	return;
340	}
341
342	gDirectory->mkdir(dir);
343	gDirectory->cd(dir);
344
345	gDirectory->mkdir("before_cuts");
346	gDirectory->mkdir("after_cuts");
347
348	hCutStatistics->Write();
349	hCutCorrelation->Write();
350
351	for (Int_t i=0; i<2; i++) {
352	if (i==0)
353	gDirectory->cd("before_cuts");
354	else
355	gDirectory->cd("after_cuts");
356
357	hNClustersITS[i] ->Write();
358	hNClustersTPC[i] ->Write();
359	hChi2PerClusterITS[i] ->Write();
360	hChi2PerClusterTPC[i] ->Write();
361
362	hC11[i] ->Write();
363	hC22[i] ->Write();
364	hC33[i] ->Write();
365	hC44[i] ->Write();
366	hC55[i] ->Write();
367
368	hDXY[i] ->Write();
369	hDZ[i] ->Write();
370	hDXYvsDZ[i] ->Write();
371
372	hDXYNormalized[i] ->Write();
373	hDZNormalized[i] ->Write();
374	hDXYvsDZNormalized[i] ->Write();
375
376	gDirectory->cd("../");
377	}
378
379	gDirectory->cd("../");
380	}
381
382
383