[u/mrichter/AliRoot.git] / PWG3 / muon / PlotEfficiency.C

/**************************************************************************
 * This file is property of and copyright by the ALICE HLT Project        *
 * All rights reserved.                                                   *
 *                                                                        *
 * Primary Authors:                                                       *
 *   Artur Szostak <artursz@iafrica.com>                                  *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/// \ingroup macros
/// \file PlotEfficiency.C
/// \brief Generates efficiency plots from the histograms generated by AliAnalysisTaskLinkToMC.
/// \author Artur Szostak <artursz@iafrica.com>
///
/// This macro is used to generate efficiency plots, fake track ratio plots
/// and also calculate the total integrated efficiency and fake track ratio.
/// The histograms generated by the AliAnalysisTaskLinkToMC analysis task is
/// used as input. The macro can be run with root as follows:
/// 
///  $ root PlotEfficiency.C\(\"hists.root\"\)
///
/// where hists.root is the name of the output file containing the histograms
/// generated by the analysis task. Note that the '\' character before the '('
/// and '"' characters is required. Alternatively run the macros as follows from
/// the root command prompt:
///
///  $ root
///  .L PlotEfficiency.C
///  PlotEfficiency("hists.root")


#if !defined(__CINT__) || defined(__MAKECINT__)
#include "Riostream.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TList.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"
#endif


void EfficiencyIntegral(const char* prefix, TH1* pass, TH1* total, double& eff, double& error, bool includeOverflow = false)
{
	/// This routine calculates the integrated efficiency for a 1D histogram.
	/// [in] \param prefix  The string to print before the calculated efficiency result.
	/// [in] \param pass  The histogram of reconstructed tracks.
	/// [in] \param total  The histogram of findable (reconstructable) tracks.
	/// [out] \param eff  Filled by the total integrated efficiency calculated.
	/// [out] \param error  Filled by the uncertainty in the efficiency value.
	/// [in] \param includeOverflow  If set then the overflow bins are also used
	///     in the calculation.
	
	TH1D passhist("temppasshist", "", 1, 0, 1);
	TH1D totalhist("temptotalhist", "", 1, 0, 1);
	if (includeOverflow)
	{
		passhist.SetBinContent(1, pass->Integral(0, pass->GetNbinsX()+1));
		totalhist.SetBinContent(1, total->Integral(0, total->GetNbinsX()+1));
	}
	else
	{
		passhist.SetBinContent(1, pass->Integral(1, pass->GetNbinsX()));
		totalhist.SetBinContent(1, total->Integral(1, total->GetNbinsX()));
	}
	TGraphAsymmErrors* effgraph = new TGraphAsymmErrors(&passhist, &totalhist);
	eff = effgraph->GetY()[0];
	double effhigh = effgraph->GetEYhigh()[0];
	double efflow = effgraph->GetEYlow()[0];
	error = effhigh > efflow ? effhigh : efflow;
	cout << prefix << eff << " + " << effhigh << " - " << efflow << endl;
}


void PlotEfficiency(const char* histfilename = "hists.root", bool includeOverflow = true)
{
	/// Opens the file containing the histograms generated by the AliAnalysisTaskLinkToMC
	/// analysis task and tries to find the first TList in the file.
	/// From the TList we try to find TH2 classes with the following names:
	/// "findableTracksHist", "foundTracksHistMC" and "fakeTracksHist"
	/// From these we plot the efficiency plots for the pT and rapidity dimension
	/// by projecting the 2D histograms. Also, the fake track ratio histograms are
	/// created. Finally the total integrated efficiency and fake track ratio is
	/// calculated and printed.
	/// \param histfilename  The name of the root file containing the histograms
	///    generated by the AliAnalysisTaskLinkToMC analysis task.
	/// \param includeOverflow  Indicates if the overflow bins should be used for
	///    for the total integral efficiency and fake track ratio calculations.
	///    The default is not to use the overflow bins.
	
	TFile file(histfilename);
	
	TList* histolist = NULL;
	TIter next(file.GetListOfKeys());
	for (TObject* key = next(); key != NULL; key = next())
	{
		TList* list = dynamic_cast<TList*>( file.Get(key->GetName()) );
		if (list != NULL)
		{
			histolist = list;
			break;
		}
	}
	
	if (histolist == NULL)
	{
		cout << "ERROR: Could not find the histogram list in file '" << histfilename << "'" << endl;
		return;
	}
	
	TH2* findableTracksHist = dynamic_cast<TH2*>( histolist->FindObject("findableTracksHist") );
	TH2* foundTracksHistMC = dynamic_cast<TH2*>( histolist->FindObject("foundTracksHistMC") );
	TH2* fakeTracksHist = dynamic_cast<TH2*>( histolist->FindObject("fakeTracksHist") );
	
	if (findableTracksHist == NULL)
	{
		cout << "ERROR: Could not find the histogram findableTracksHist in TList '" << histolist->GetName() << "'" << endl;
		return;
	}
	if (foundTracksHistMC == NULL)
	{
		cout << "ERROR: Could not find the histogram foundTracksHistMC in TList '" << histolist->GetName() << "'" << endl;
		return;
	}
	if (fakeTracksHist == NULL)
	{
		cout << "ERROR: Could not find the histogram fakeTracksHist in TList '" << histolist->GetName() << "'" << endl;
		return;
	}
	
	TH1* projeffx = foundTracksHistMC->ProjectionX();
	TH1* projfindx = findableTracksHist->ProjectionX();
	TGraphAsymmErrors* effpt = new TGraphAsymmErrors(projeffx, projfindx);
	effpt->Draw("ap");
	effpt->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
	effpt->GetHistogram()->SetXTitle(projeffx->GetXaxis()->GetTitle());
	effpt->GetHistogram()->SetYTitle("#epsilon");
	effpt->GetYaxis()->SetRangeUser(0, 1.1);
	
	new TCanvas;
	TH1* projfakex = fakeTracksHist->ProjectionX();
	TGraphAsymmErrors* fakept = new TGraphAsymmErrors(projfakex, projfindx);
	fakept->Draw("ap");
	fakept->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
	fakept->GetHistogram()->SetXTitle(projfakex->GetXaxis()->GetTitle());
	fakept->GetHistogram()->SetYTitle("fake tracks / findable tracks");
	
	new TCanvas;
	TH1* projeffy = foundTracksHistMC->ProjectionY();
	TH1* projfindy = findableTracksHist->ProjectionY();
	TGraphAsymmErrors* effy = new TGraphAsymmErrors(projeffy, projfindy);
	effy->Draw("ap");
	effy->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
	effy->GetHistogram()->SetXTitle(projeffy->GetXaxis()->GetTitle());
	effy->GetHistogram()->SetYTitle("#epsilon");
	effy->GetYaxis()->SetRangeUser(0, 1.1);
	
	new TCanvas;
	TH1* projfakey = fakeTracksHist->ProjectionY();
	TGraphAsymmErrors* fakey = new TGraphAsymmErrors(projfakey, projfindy);
	fakey->Draw("ap");
	fakey->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
	fakey->GetHistogram()->SetXTitle(projfakey->GetXaxis()->GetTitle());
	fakey->GetHistogram()->SetYTitle("fake tracks / findable tracks");
	
	if (includeOverflow)
	{
		cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
		cout << "found tracks    = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
		cout << "fake tracks     = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
	}
	else
	{
		cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
		cout << "found tracks    = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
		cout << "fake tracks     = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
	}
	double eff, fake, deff, dfake;
	EfficiencyIntegral("Total efficiency = ", projeffx, projfindx, eff, deff, includeOverflow);
	EfficiencyIntegral("Total fake ratio = ", projfakex, projfindx, fake, dfake, includeOverflow);
	cout << "Sum = " << eff + fake << " +/- " << sqrt(deff*deff + dfake*dfake) << endl;
}
Commit	Line	Data
4c904d93	1	/**************************************************************************
	2	* This file is property of and copyright by the ALICE HLT Project *
	3	* All rights reserved. *
	4	* *
	5	* Primary Authors: *
	6	* Artur Szostak <artursz@iafrica.com> *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16
	17	/// \ingroup macros
	18	/// \file PlotEfficiency.C
	19	/// \brief Generates efficiency plots from the histograms generated by AliAnalysisTaskLinkToMC.
	20	/// \author Artur Szostak <artursz@iafrica.com>
	21	///
	22	/// This macro is used to generate efficiency plots, fake track ratio plots
	23	/// and also calculate the total integrated efficiency and fake track ratio.
	24	/// The histograms generated by the AliAnalysisTaskLinkToMC analysis task is
	25	/// used as input. The macro can be run with root as follows:
	26	///
	27	/// $ root PlotEfficiency.C\(\"hists.root\"\)
	28	///
	29	/// where hists.root is the name of the output file containing the histograms
	30	/// generated by the analysis task. Note that the '\' character before the '('
	31	/// and '"' characters is required. Alternatively run the macros as follows from
	32	/// the root command prompt:
	33	///
	34	/// $ root
	35	/// .L PlotEfficiency.C
	36	/// PlotEfficiency("hists.root")
	37
	38
	39	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	40	#include "Riostream.h"
	41	#include "TH1.h"
	42	#include "TH2.h"
	43	#include "TFile.h"
	44	#include "TList.h"
	45	#include "TGraphAsymmErrors.h"
	46	#include "TCanvas.h"
	47	#endif
	48
	49
	50	void EfficiencyIntegral(const char* prefix, TH1* pass, TH1* total, double& eff, double& error, bool includeOverflow = false)
	51	{
	52	/// This routine calculates the integrated efficiency for a 1D histogram.
	53	/// [in] \param prefix The string to print before the calculated efficiency result.
	54	/// [in] \param pass The histogram of reconstructed tracks.
	55	/// [in] \param total The histogram of findable (reconstructable) tracks.
	56	/// [out] \param eff Filled by the total integrated efficiency calculated.
	57	/// [out] \param error Filled by the uncertainty in the efficiency value.
	58	/// [in] \param includeOverflow If set then the overflow bins are also used
	59	/// in the calculation.
	60
	61	TH1D passhist("temppasshist", "", 1, 0, 1);
	62	TH1D totalhist("temptotalhist", "", 1, 0, 1);
	63	if (includeOverflow)
	64	{
65	passhist.SetBinContent(1, pass->Integral(0, pass->GetNbinsX()+1));
66	totalhist.SetBinContent(1, total->Integral(0, total->GetNbinsX()+1));
67	}
68	else
69	{
70	passhist.SetBinContent(1, pass->Integral(1, pass->GetNbinsX()));
71	totalhist.SetBinContent(1, total->Integral(1, total->GetNbinsX()));
72	}
73	TGraphAsymmErrors* effgraph = new TGraphAsymmErrors(&passhist, &totalhist);
74	eff = effgraph->GetY()[0];
75	double effhigh = effgraph->GetEYhigh()[0];
76	double efflow = effgraph->GetEYlow()[0];
77	error = effhigh > efflow ? effhigh : efflow;
78	cout << prefix << eff << " + " << effhigh << " - " << efflow << endl;
79	}
80
81
82	void PlotEfficiency(const char* histfilename = "hists.root", bool includeOverflow = true)
83	{
84	/// Opens the file containing the histograms generated by the AliAnalysisTaskLinkToMC
85	/// analysis task and tries to find the first TList in the file.
86	/// From the TList we try to find TH2 classes with the following names:
87	/// "findableTracksHist", "foundTracksHistMC" and "fakeTracksHist"
88	/// From these we plot the efficiency plots for the pT and rapidity dimension
89	/// by projecting the 2D histograms. Also, the fake track ratio histograms are
90	/// created. Finally the total integrated efficiency and fake track ratio is
91	/// calculated and printed.
92	/// \param histfilename The name of the root file containing the histograms
93	/// generated by the AliAnalysisTaskLinkToMC analysis task.
94	/// \param includeOverflow Indicates if the overflow bins should be used for
95	/// for the total integral efficiency and fake track ratio calculations.
96	/// The default is not to use the overflow bins.
97
98	TFile file(histfilename);
99
100	TList* histolist = NULL;
101	TIter next(file.GetListOfKeys());
102	for (TObject* key = next(); key != NULL; key = next())
103	{
104	TList* list = dynamic_cast<TList*>( file.Get(key->GetName()) );
105	if (list != NULL)
106	{
107	histolist = list;
108	break;
109	}
110	}
111
112	if (histolist == NULL)
113	{
114	cout << "ERROR: Could not find the histogram list in file '" << histfilename << "'" << endl;
115	return;
116	}
117
118	TH2* findableTracksHist = dynamic_cast<TH2*>( histolist->FindObject("findableTracksHist") );
119	TH2* foundTracksHistMC = dynamic_cast<TH2*>( histolist->FindObject("foundTracksHistMC") );
120	TH2* fakeTracksHist = dynamic_cast<TH2*>( histolist->FindObject("fakeTracksHist") );
121
122	if (findableTracksHist == NULL)
123	{
124	cout << "ERROR: Could not find the histogram findableTracksHist in TList '" << histolist->GetName() << "'" << endl;
125	return;
126	}
127	if (foundTracksHistMC == NULL)
128	{
129	cout << "ERROR: Could not find the histogram foundTracksHistMC in TList '" << histolist->GetName() << "'" << endl;
130	return;
131	}
132	if (fakeTracksHist == NULL)
133	{
134	cout << "ERROR: Could not find the histogram fakeTracksHist in TList '" << histolist->GetName() << "'" << endl;
135	return;
136	}
137
138	TH1* projeffx = foundTracksHistMC->ProjectionX();
139	TH1* projfindx = findableTracksHist->ProjectionX();
140	TGraphAsymmErrors* effpt = new TGraphAsymmErrors(projeffx, projfindx);
141	effpt->Draw("ap");
142	effpt->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
143	effpt->GetHistogram()->SetXTitle(projeffx->GetXaxis()->GetTitle());
144	effpt->GetHistogram()->SetYTitle("#epsilon");
145	effpt->GetYaxis()->SetRangeUser(0, 1.1);
146
147	new TCanvas;
148	TH1* projfakex = fakeTracksHist->ProjectionX();
149	TGraphAsymmErrors* fakept = new TGraphAsymmErrors(projfakex, projfindx);
150	fakept->Draw("ap");
151	fakept->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
152	fakept->GetHistogram()->SetXTitle(projfakex->GetXaxis()->GetTitle());
153	fakept->GetHistogram()->SetYTitle("fake tracks / findable tracks");
154
155	new TCanvas;
156	TH1* projeffy = foundTracksHistMC->ProjectionY();
157	TH1* projfindy = findableTracksHist->ProjectionY();
158	TGraphAsymmErrors* effy = new TGraphAsymmErrors(projeffy, projfindy);
159	effy->Draw("ap");
160	effy->GetHistogram()->SetTitle("Efficiency of track reconstruction.");
161	effy->GetHistogram()->SetXTitle(projeffy->GetXaxis()->GetTitle());
162	effy->GetHistogram()->SetYTitle("#epsilon");
163	effy->GetYaxis()->SetRangeUser(0, 1.1);
164
165	new TCanvas;
166	TH1* projfakey = fakeTracksHist->ProjectionY();
167	TGraphAsymmErrors* fakey = new TGraphAsymmErrors(projfakey, projfindy);
168	fakey->Draw("ap");
169	fakey->GetHistogram()->SetTitle("Fake reconstructed tracks ratio.");
170	fakey->GetHistogram()->SetXTitle(projfakey->GetXaxis()->GetTitle());
171	fakey->GetHistogram()->SetYTitle("fake tracks / findable tracks");
172
173	if (includeOverflow)
174	{
175	cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
176	cout << "found tracks = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
177	cout << "fake tracks = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
178	}
179	else
180	{
181	cout << "findable tracks = " << projfindx->Integral(0, projfindx->GetNbinsX()+1) << endl;
182	cout << "found tracks = " << projeffx->Integral(0, projeffx->GetNbinsX()+1) << endl;
183	cout << "fake tracks = " << projfakex->Integral(0, projfakex->GetNbinsX()+1) << endl;
184	}
185	double eff, fake, deff, dfake;
186	EfficiencyIntegral("Total efficiency = ", projeffx, projfindx, eff, deff, includeOverflow);
187	EfficiencyIntegral("Total fake ratio = ", projfakex, projfindx, fake, dfake, includeOverflow);
188	cout << "Sum = " << eff + fake << " +/- " << sqrt(deffdeff + dfakedfake) << endl;
189	}
190