#include "TROOT.h"
#include "TChain.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TString.h"
#include "TLatex.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TLegend.h"

#include <vector>
#include <list>
#include <string>

TCanvas* canvas =  new TCanvas("can","can",1000,1000);

int findXvalue(TH1D* projection, double splitFraction, bool signFlip){

	if(signFlip)splitFraction = 1. - splitFraction;


	if( splitFraction <=0 || splitFraction >=1 ){
		std::cout << "ERROR -  SplitFraction: " << splitFraction << std::endl; 
		return -1;
	}

	double*  integral = projection->GetIntegral();
	int i(0);
	while( integral[i]< splitFraction ){
		++i;
	}
	return i;
}

double findFraction(TH1D* projection, int bin, bool signFlip){

	if (bin < 0 || bin  > projection->GetNbinsX() ){
		std::cout << "ERROR -  BinCrazy: " << bin << std::endl; 
		return 0;
	}

	double*  integral = projection->GetIntegral();
	double result = integral[bin];
	if(signFlip){
		result = 1 - result;
	}

	return result;
}

void printGraphInclusiveFakeRate(TH3D* histo3dNpartX, TH3D* histo3dNpart2, double splitFraction) {
	
	int nToTpoints = histo3dNpartX->GetYaxis()->GetNbins();
	int nLayers = histo3dNpartX->GetXaxis()->GetNbins();
	cout << "nLayers = " << nLayers << endl;
	std::vector<TGraphErrors*> graphs(nLayers);
	TLegend* l=new TLegend(0.11, 0.89, 0.26, 0.68);
	
	for (int layer(1); layer<=nLayers; ++layer) {
		TGraphErrors *g = new TGraphErrors(histo3dNpartX->GetYaxis()->GetNbins());
		TString gname("Layer ");
		if (layer==1) gname += "IBL";
		else gname+=Form("%d", layer-1);
		g->SetName(gname);
		
		for (int ToTp(1); ToTp<=nToTpoints; ++ToTp) {
		
			TH1D* histo1 = histo3dNpartX->ProjectionZ("histo1", layer, layer, ToTp, ToTp);
			TH1D* histo2 = histo3dNpart2->ProjectionZ("histo2", layer, layer, ToTp, ToTp);
	
			double fakerate(0.);
			int cutInt(0);
	
			if (histo1->GetEntries() < 100 || histo2->GetEntries() < 100) {
				cout << "Nr of entries too small, setting fake rate to zero" << endl;
				// cutInt = 0;
				fakerate = 0.;
			}
		
			else {
				cutInt  =  findXvalue(  histo2 , splitFraction, false );

				if( cutInt > 0){
				// with efficiency = splitFraction, frac gives the fraction of 1-part splits
					fakerate =  findFraction(  histo1 , cutInt, false );
				}
			}
			if (fakerate>0.) g->SetPoint(ToTp-1, histo3dNpartX->GetYaxis()->GetBinCenter(ToTp), fakerate);
			delete histo1;
			delete histo2;
		}
		l->AddEntry(g, gname, "p");
// 		graphs[layer-1]=g;
		graphs[layer-1]=g;
	}
	for (unsigned int i(0); i<graphs.size(); ++i) {
		if (i==0) graphs[i]->Draw("AP");
		else graphs[i]->Draw("P");
		graphs[i]->SetMarkerStyle(8);
		graphs[i]->SetMarkerColor(i+1);
		graphs[i]->SetMarkerSize(1.1);
		
	}
	l->Draw();
	canvas->Print("FakeRateInclusiveVsToTtreshmod.png");
	
	for (unsigned int i(0); i<graphs.size(); ++i) {
		delete graphs[i];//[i];
	}	
}

void makeplot(TString fname="perf.root")
{

	gStyle->SetOptStat(0);
	
	//Open file
	TFile* file =  new TFile(fname);
	
	//Get 2 histograms
	TH3D* SplitProb1vsLayerToTtreshmodNpart1 = (TH3D*) file->Get("SplitProb1vsLayerToTtreshmodNpart1");
	TH3D* SplitProb1vsLayerToTtreshmodNpart2 = (TH3D*) file->Get("SplitProb1vsLayerToTtreshmodNpart2");
	
	//Parameter needed for plotting
	double splitFraction = 0.75;
	//Pass these histograms to a function which processes their content and plots some info
	printGraphInclusiveFakeRate(SplitProb1vsLayerToTtreshmodNpart1, SplitProb1vsLayerToTtreshmodNpart2, splitFraction);
	
}