/**************************************************************************
 * 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;
}