[u/mrichter/AliRoot.git] / PWGLF / totEt / macros / emEt / EMCal / CalcEffic.C

/*
 *  CalcEffic.C
 */

TEfficiency* CalcEffic(const char* fileName, Bool_t makeDraw=kFALSE)
{
	// open the input file
	TFile *file = new TFile(fileName);
	
	// ********************************************
	// parameters to check before running the macro
	// ********************************************

	const Float_t E_MIN = 0.1;
	const Int_t NHISTS = 6; // Check the number of histograms for different particle species
	const Int_t NOUTPUTS = 3;
	const Int_t NHISTOUT[NOUTPUTS] = {6,3,3};
	Int_t IHISTOUT[NOUTPUTS][NHISTS] = {{0,1,2,3,4,5},{0,1,2,-1,-1,-1},{3,4,5,-1,-1,-1}};
	
	Int_t style[NOUTPUTS] = {20,21,22};
	Int_t color[NOUTPUTS] = {1,2,4};
		
	const Int_t fgNumOfEBins = 78; // Check the number of eta bins in the histograms
	const Int_t fgNumOfEtaBins = 16; // Check the number of E bins in the histograms
	Double_t fgEAxis[79]={0., 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
		1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8.,8.5,9.,9.5,10.,11.,
		12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29.,30.,31.,
		32.,33.,34.,35.,36.,37.,38.,39.,40.,41.,42.,43.,44.,45.,46.,47.,48.,49.,50.};
	
	// declare histograms and graphs
	TH2F *histNum[NHISTS];
	TH2F *histDen[NHISTS];
	TGraphErrors *graph[NOUTPUTS];
	TGraphErrors *graphNum[NOUTPUTS];
	TGraphErrors *graphDen[NOUTPUTS];
	TEfficiency *effic[NOUTPUTS];
	char efficName[50];
	
	//define canvas
	gStyle->SetOptTitle(0);
	gStyle->SetOptStat(0);
	gStyle->SetOptFit(0);	
	
	TCanvas *c1;
	TCanvas *c2;
	
	if (makeDraw)
	{
		c1 = new TCanvas("c1","c1",500,400);
		//c1->SetTopMargin(0.04);
		//c1->SetRightMargin(0.04);
		//c1->SetLeftMargin(0.181452);
		//c1->SetBottomMargin(0.134409);
		c1->SetBorderSize(0);
		c1->SetFillColor(0);
		c1->SetBorderMode(0);
		c1->SetFrameFillColor(0);
		c1->SetFrameBorderMode(0);
		
		c2 = new TCanvas("c2","c2",500,400);
		c2->Divide(1,3);
	}
	
	// retrieve the input list of histogram. Check the TList name in the input file.
	TList *list = (TList*) file->Get("out1");
	
	// retrieve the histograms in the list. Check the name of the histograms
	histNum[0] = (TH2F*)list->FindObject("fHistElectronRec_EtaE_EmcalMC");
	histNum[1] = (TH2F*)list->FindObject("fHistConvElectronRec_EtaE_EmcalMC");
	histNum[2] = (TH2F*)list->FindObject("fHistScatElectronRec_EtaE_EmcalMC");
	histNum[3] = (TH2F*)list->FindObject("fHistGammaRec_EtaE_EmcalMC");
	histNum[4] = (TH2F*)list->FindObject("fHistAnnihGammaRec_EtaE_EmcalMC");
	histNum[5] = (TH2F*)list->FindObject("fHistScatGammaRec_EtaE_EmcalMC");
	
	// retrieve the histograms in the list. Check the name of the histograms
	histDen[0] = (TH2F*)list->FindObject("fHistElectronAcc_EtaE_EmcalMC");
	histDen[1] = (TH2F*)list->FindObject("fHistConvElectronAcc_EtaE_EmcalMC");
	histDen[2] = (TH2F*)list->FindObject("fHistScatElectronAcc_EtaE_EmcalMC");
	histDen[3] = (TH2F*)list->FindObject("fHistGammaAcc_EtaE_EmcalMC");
	histDen[4] = (TH2F*)list->FindObject("fHistAnnihGammaAcc_EtaE_EmcalMC");
	histDen[5] = (TH2F*)list->FindObject("fHistScatGammaAcc_EtaE_EmcalMC");
	
	// ********************************************

	Float_t x[fgNumOfEBins]={0}, ex[fgNumOfEBins]={0};
	Float_t y[fgNumOfEBins]={0}, ey[fgNumOfEBins]={0};
	Float_t num[fgNumOfEBins]={0}, den[fgNumOfEBins]={0};
	
	// loop over different desired outputs
	for (int iOut=0; iOut<NOUTPUTS; iOut++)
	{
		sprintf(efficName,"effic_%d",iOut);
		effic[iOut] = new TEfficiency(efficName,efficName,fgNumOfEBins,fgEAxis);

		// loop over E bins
		for (int ix=0; ix<fgNumOfEBins; ix++)
		{
			//check minimum energy
			if (histNum[0]->GetXaxis()->GetBinLowEdge(ix+1) < E_MIN)
				continue;
			
			// initialize ET variables for a new particle species
			x[ix]=histNum[0]->GetXaxis()->GetBinCenter(ix+1);
			y[ix]=0;
			ex[ix]=0;
			ey[ix]=0;
			num[ix] = 0;
			den[ix] = 0;
			
			// loop over eta bins
			for (int iy=0; iy<fgNumOfEtaBins; iy++)
			{
				for (int iHist=0; iHist<NHISTOUT[iOut]; iHist++)
				{
					num[ix] += histNum[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
					den[ix] += histDen[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
				}
			}
			
			if ((num[ix]>0) && (den[ix]>0))
			{
				effic[iOut]->SetTotalEvents(ix,den[ix]);
				effic[iOut]->SetPassedEvents(ix,num[ix]);

				y[ix] = num[ix]/den[ix];
				ey[ix] = y[ix]*sqrt(1/num[ix]+1/den[ix]);				
			}
			else
			{
				y[ix] = 0;	
				ey[ix] = 0;
			}
			
		} // end of loop over E bins

		graph[iOut] = new TGraphErrors(fgNumOfEBins,x,y,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut
		graphNum[iOut] = new TGraphErrors(fgNumOfEBins,x,num,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut	
		graphDen[iOut] = new TGraphErrors(fgNumOfEBins,x,den,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut

	} // end of loop over different outputs

	
	// Draw the plot
	
	if (makeDraw)
	{
		for (int i=0; i<NOUTPUTS; i++)
		{		
			c2->cd(i);
			
			graphDen[i]->SetMarkerStyle(style[i]);
			graphDen[i]->SetMarkerColor(color[i]);
			graphDen[i]->SetLineColor(color[i]);
			graphDen[i]->SetFillColor(0);
			if (i == 0) 
			{
				graphDen[i]->GetXaxis()->SetTitle("E (GeV)");
				graphDen[i]->GetYaxis()->SetTitle("effic");
				//graphDen[i]->SetMaximum(1.0);
				graphDen[i]->SetMinimum(0.0);
				graphDen[i]->Draw("AP");
			}
			else 
				graphDen[i]->Draw("P");
			
			graphNum[i]->SetMarkerStyle(style[i]+4);
			graphNum[i]->SetMarkerColor(color[i]);
			graphNum[i]->SetLineColor(color[i]);
			graphNum[i]->SetFillColor(0);
			graphNum[i]->Draw("P");		
		}
		
		c1->cd();
		/*
		for (int i=0; i<NOUTPUTS; i++)
		{
			effic[i]->SetMarkerStyle(style[i]);
			effic[i]->SetMarkerColor(color[i]);
			effic[i]->SetLineColor(color[i]);
			effic[i]->SetFillColor(0);
			effic[i]->SetTitle("efficiency;E (GeV); #epsilon");
			if (i == 0) 
			{
				effic[i]->Draw("AP");
			}
			else 
				effic[i]->Draw("Psame");
		}
		*/
		for (int i=0; i<NOUTPUTS; i++)
		{
			graph[i]->SetMarkerStyle(style[i]);
			graph[i]->SetMarkerColor(color[i]);
			graph[i]->SetLineColor(color[i]);
			graph[i]->SetFillColor(0);
			if (i == 0) 
			{
				graph[i]->Draw("AP");
			}
			else 
				graph[i]->Draw("P");
		}
		
		TLegend *leg = new TLegend(0.65,0.2,0.95,0.5);
		leg->AddEntry(effic[0],"electrons+gammas");
		leg->AddEntry(effic[1],"electrons");
		leg->AddEntry(effic[2],"gammas");
		leg->SetFillStyle(0);
		leg->SetFillColor(0);
		leg->SetBorderSize(0);
		leg->SetTextSize(0.03);
		leg->Draw();
	}
	
	return effic[0];
}
Commit	Line	Data
f18f3638	1	/*
	2	* CalcEffic.C
	3	*/
	4
	5	TEfficiency* CalcEffic(const char* fileName, Bool_t makeDraw=kFALSE)
	6	{
	7	// open the input file
	8	TFile *file = new TFile(fileName);
	9
	10	// ********************************************
	11	// parameters to check before running the macro
	12	// ********************************************
	13
	14	const Float_t E_MIN = 0.1;
	15	const Int_t NHISTS = 6; // Check the number of histograms for different particle species
	16	const Int_t NOUTPUTS = 3;
	17	const Int_t NHISTOUT[NOUTPUTS] = {6,3,3};
	18	Int_t IHISTOUT[NOUTPUTS][NHISTS] = {{0,1,2,3,4,5},{0,1,2,-1,-1,-1},{3,4,5,-1,-1,-1}};
	19
	20	Int_t style[NOUTPUTS] = {20,21,22};
	21	Int_t color[NOUTPUTS] = {1,2,4};
	22
	23	const Int_t fgNumOfEBins = 78; // Check the number of eta bins in the histograms
	24	const Int_t fgNumOfEtaBins = 16; // Check the number of E bins in the histograms
	25	Double_t fgEAxis[79]={0., 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
	26	1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,6.5,7.,7.5,8.,8.5,9.,9.5,10.,11.,
	27	12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29.,30.,31.,
	28	32.,33.,34.,35.,36.,37.,38.,39.,40.,41.,42.,43.,44.,45.,46.,47.,48.,49.,50.};
	29
	30	// declare histograms and graphs
	31	TH2F *histNum[NHISTS];
	32	TH2F *histDen[NHISTS];
	33	TGraphErrors *graph[NOUTPUTS];
	34	TGraphErrors *graphNum[NOUTPUTS];
	35	TGraphErrors *graphDen[NOUTPUTS];
	36	TEfficiency *effic[NOUTPUTS];
	37	char efficName[50];
	38
	39	//define canvas
	40	gStyle->SetOptTitle(0);
	41	gStyle->SetOptStat(0);
	42	gStyle->SetOptFit(0);
	43
	44	TCanvas *c1;
	45	TCanvas *c2;
	46
	47	if (makeDraw)
	48	{
	49	c1 = new TCanvas("c1","c1",500,400);
	50	//c1->SetTopMargin(0.04);
	51	//c1->SetRightMargin(0.04);
	52	//c1->SetLeftMargin(0.181452);
	53	//c1->SetBottomMargin(0.134409);
	54	c1->SetBorderSize(0);
	55	c1->SetFillColor(0);
	56	c1->SetBorderMode(0);
	57	c1->SetFrameFillColor(0);
	58	c1->SetFrameBorderMode(0);
	59
	60	c2 = new TCanvas("c2","c2",500,400);
	61	c2->Divide(1,3);
	62	}
	63
	64	// retrieve the input list of histogram. Check the TList name in the input file.
65	TList list = (TList) file->Get("out1");
66
67	// retrieve the histograms in the list. Check the name of the histograms
68	histNum[0] = (TH2F*)list->FindObject("fHistElectronRec_EtaE_EmcalMC");
69	histNum[1] = (TH2F*)list->FindObject("fHistConvElectronRec_EtaE_EmcalMC");
70	histNum[2] = (TH2F*)list->FindObject("fHistScatElectronRec_EtaE_EmcalMC");
71	histNum[3] = (TH2F*)list->FindObject("fHistGammaRec_EtaE_EmcalMC");
72	histNum[4] = (TH2F*)list->FindObject("fHistAnnihGammaRec_EtaE_EmcalMC");
73	histNum[5] = (TH2F*)list->FindObject("fHistScatGammaRec_EtaE_EmcalMC");
74
75	// retrieve the histograms in the list. Check the name of the histograms
76	histDen[0] = (TH2F*)list->FindObject("fHistElectronAcc_EtaE_EmcalMC");
77	histDen[1] = (TH2F*)list->FindObject("fHistConvElectronAcc_EtaE_EmcalMC");
78	histDen[2] = (TH2F*)list->FindObject("fHistScatElectronAcc_EtaE_EmcalMC");
79	histDen[3] = (TH2F*)list->FindObject("fHistGammaAcc_EtaE_EmcalMC");
80	histDen[4] = (TH2F*)list->FindObject("fHistAnnihGammaAcc_EtaE_EmcalMC");
81	histDen[5] = (TH2F*)list->FindObject("fHistScatGammaAcc_EtaE_EmcalMC");
82
83	// ********************************************
84
85	Float_t x[fgNumOfEBins]={0}, ex[fgNumOfEBins]={0};
86	Float_t y[fgNumOfEBins]={0}, ey[fgNumOfEBins]={0};
87	Float_t num[fgNumOfEBins]={0}, den[fgNumOfEBins]={0};
88
89	// loop over different desired outputs
90	for (int iOut=0; iOut<NOUTPUTS; iOut++)
91	{
92	sprintf(efficName,"effic_%d",iOut);
93	effic[iOut] = new TEfficiency(efficName,efficName,fgNumOfEBins,fgEAxis);
94
95	// loop over E bins
96	for (int ix=0; ix<fgNumOfEBins; ix++)
97	{
98	//check minimum energy
99	if (histNum[0]->GetXaxis()->GetBinLowEdge(ix+1) < E_MIN)
100	continue;
101
102	// initialize ET variables for a new particle species
103	x[ix]=histNum[0]->GetXaxis()->GetBinCenter(ix+1);
104	y[ix]=0;
105	ex[ix]=0;
106	ey[ix]=0;
107	num[ix] = 0;
108	den[ix] = 0;
109
110	// loop over eta bins
111	for (int iy=0; iy<fgNumOfEtaBins; iy++)
112	{
113	for (int iHist=0; iHist<NHISTOUT[iOut]; iHist++)
114	{
115	num[ix] += histNum[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
116	den[ix] += histDen[IHISTOUT[iOut][iHist]]->GetBinContent(ix+1,iy+1); // sum over all E bins in order to get total ET
117	}
118	}
119
120	if ((num[ix]>0) && (den[ix]>0))
121	{
122	effic[iOut]->SetTotalEvents(ix,den[ix]);
123	effic[iOut]->SetPassedEvents(ix,num[ix]);
124
125	y[ix] = num[ix]/den[ix];
126	ey[ix] = y[ix]*sqrt(1/num[ix]+1/den[ix]);
127	}
128	else
129	{
130	y[ix] = 0;
131	ey[ix] = 0;
132	}
133
134	} // end of loop over E bins
135
136	graph[iOut] = new TGraphErrors(fgNumOfEBins,x,y,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut
137	graphNum[iOut] = new TGraphErrors(fgNumOfEBins,x,num,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut
138	graphDen[iOut] = new TGraphErrors(fgNumOfEBins,x,den,ex,ey); // graphic of ET(>E_cut)/ET(total) for a given particle species and E cut
139
140	} // end of loop over different outputs
141
142
143	// Draw the plot
144
145	if (makeDraw)
146	{
147	for (int i=0; i<NOUTPUTS; i++)
148	{
149	c2->cd(i);
150
151	graphDen[i]->SetMarkerStyle(style[i]);
152	graphDen[i]->SetMarkerColor(color[i]);
153	graphDen[i]->SetLineColor(color[i]);
154	graphDen[i]->SetFillColor(0);
155	if (i == 0)
156	{
157	graphDen[i]->GetXaxis()->SetTitle("E (GeV)");
158	graphDen[i]->GetYaxis()->SetTitle("effic");
159	//graphDen[i]->SetMaximum(1.0);
160	graphDen[i]->SetMinimum(0.0);
161	graphDen[i]->Draw("AP");
162	}
163	else
164	graphDen[i]->Draw("P");
165
166	graphNum[i]->SetMarkerStyle(style[i]+4);
167	graphNum[i]->SetMarkerColor(color[i]);
168	graphNum[i]->SetLineColor(color[i]);
169	graphNum[i]->SetFillColor(0);
170	graphNum[i]->Draw("P");
171	}
172
173	c1->cd();
174	/*
175	for (int i=0; i<NOUTPUTS; i++)
176	{
177	effic[i]->SetMarkerStyle(style[i]);
178	effic[i]->SetMarkerColor(color[i]);
179	effic[i]->SetLineColor(color[i]);
180	effic[i]->SetFillColor(0);
181	effic[i]->SetTitle("efficiency;E (GeV); #epsilon");
182	if (i == 0)
183	{
184	effic[i]->Draw("AP");
185	}
186	else
187	effic[i]->Draw("Psame");
188	}
189	*/
190	for (int i=0; i<NOUTPUTS; i++)
191	{
192	graph[i]->SetMarkerStyle(style[i]);
193	graph[i]->SetMarkerColor(color[i]);
194	graph[i]->SetLineColor(color[i]);
195	graph[i]->SetFillColor(0);
196	if (i == 0)
197	{
198	graph[i]->Draw("AP");
199	}
200	else
201	graph[i]->Draw("P");
202	}
203
204	TLegend *leg = new TLegend(0.65,0.2,0.95,0.5);
205	leg->AddEntry(effic[0],"electrons+gammas");
206	leg->AddEntry(effic[1],"electrons");
207	leg->AddEntry(effic[2],"gammas");
208	leg->SetFillStyle(0);
209	leg->SetFillColor(0);
210	leg->SetBorderSize(0);
211	leg->SetTextSize(0.03);
212	leg->Draw();
213	}
214
215	return effic[0];
216	}