#include<iostream>
#include<vector>
#include<string>
#include "TString.h"
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"

int BeamProfile(TString RootFile,TString RecoFile, Int_t name)
{
//	gROOT->SetStyle("Plain");
//    	gStyle->SetOptFit(1111);
//        gStyle->SetOptStat(1111);
//    	gStyle->SetPalette(1);        

	gROOT->ProcessLine(".L tdrstyle.C");
	gROOT->ProcessLine(".x rootlogon.C");
	setTDRStyle();

	TGaxis::SetMaxDigits(3);

	cmsprem = new TLatex(0,101,"CMS Preliminary");
	cmsprem->SetTextSize(0.04);

	  const Int_t NRGBs = 5;
	  const Int_t NCont = 255;
	  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
	  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
	  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
	  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
	  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
      	  gStyle->SetNumberContours(NCont);

	std::vector<TFile*> InputFiles;
	InputFiles.clear();
	
	InputFiles.push_back(TFile::Open(RootFile));

	const Int_t NumOfFiles = InputFiles.size();
	cout<<"Num of FIles = "<<NumOfFiles<<endl;

	for(unsigned int j = 0; j < InputFiles.size(); ++j ) {
	TTree * tmpTree = (TTree*)InputFiles.at(j)->Get("rd51tbgeo");
	tmpTree->AddFriend("trackertree",RecoFile);       

	TCanvas *canvas_prof = new TCanvas("canvas_prof","canvas_prof",1);
	gStyle->SetCanvasDefH(400);
	gStyle->SetCanvasDefW(700);
	//gPad->SetLeftMargin(0.15); 
	//gPad->SetRightMargin(-0.54); 
//	canvas_prof->Divide(2,2);
	cout<<"###################"<<endl;
//	canvas_prof->cd(1);
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
//		CORRELATION PLOTS
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
const Int_t nbranch = 12;
char *bnames[nbranch]={	"g1xcl.geoposX:g1xcl.geoposX",	"g1xcl.geoposX:g2xcl.geoposX",	"g1xcl.geoposX:g3xcl.geoposX",
			"g1ycl.geoposY:g1ycl.geoposY",	"g1ycl.geoposY:g2ycl.geoposY",	"g1ycl.geoposY:g3ycl.geoposY",
			"g1ycl.geoposY:sCMSNS2LC1.geoposX",	"g1ycl.geoposY:sCMSNS2LC2.geoposX",
			"g1xcl.geoposX:sCMSNS2LC1.geoposX",	"g1xcl.geoposX:sCMSNS2LC2.geoposX",
			"g1ycl.geoposY:sCMSNS2LC3.geoposX",	"g1xcl.geoposX:sCMSNS2LC3.geoposX"};

char *fnames[nbranch] = {"g1x_geoposX_vs_g1x_geoposX",	"g1x_geoposX_vs_g2x_geoposX",	"g1x_geoposX_vs_g3x_geoposX",
			"g1y_geoposY_vs_g1y_geoposY",	"g1y_geoposY_vs_g2y_geoposY",	"g1y_geoposY_vs_g3y_geoposY",
			"g1y_geoposY_vs_sCMSLC1_geoposX",	"g1y_geoposY_vs_sCMSLC2_geoposX",
			"g1x_geoposX_vs_sCMSLC1_geoposX",	"g1x_geoposX_vs_sCMSLC2_geoposX",
			"g1y_geoposY_vs_sCMSLC3_geoposX",	"g1x_geoposX_vs_sCMSLC3_geoposX"};
const Float_t range[4*nbranch] ={0,100,0,100,		-10,110,-10,110,		-10,110,-10,110,
				0,100,0,100,		-10,110,-10,110,		-10,110,-10,110,
				-10,110,-10,110,		-10,110,-10,110,
				-10,110,-10,110,		-10,110,-10,110,
				-10,110,-10,110,		-10,110,-10,110};

char *title[2*nbranch] = {"g1x [mm]","g1x [mm]",		"g1x [mm]","g2x [mm]",	"g1x [mm]","g3x [mm]",
			   "g1y [mm]","g1y [mm]",		"g1y [mm]","g2y [mm]",	"g1y [mm]","g3y [mm]",
			   "g1y [mm]","GE11_1(x) [mm]",	"g1y [mm]","GE11_2(x) [mm]",	
			   "g1x [mm]","GE11_1(x) [mm]",	"g1x [mm]","GE11_2(x) [mm]",	
			   "g1y [mm]","GE11_3(x) [mm]",	"g1x [mm]","GE11_3(x) [mm]"};


TH2F * hist[nbranch];
TProfile * hprofile[nbranch];
TCanvas *Canvas[nbranch];
TF1 *function[nbranch];

Char_t message[80];
vector<int> vecm,vecn; 
Double_t fit_Low, fit_High;
Int_t k=0;

ofstream o_file;
o_file.open("fit_detail.txt");

o_file<<"NAME_OF_PROFILE_PLOT\t\t\tChi-sqr/NDF\tINTERCEPT(P0)\tINTERCEPT_ERROR\t\tSLOPE(P1)\tSLOPE_ERROR"<<endl;
for(Int_t i=0;i<nbranch;i++){
	k=i*4;	m=i*2;
	Canvas[i] = new TCanvas(Form("Canvas%d",i));
	Canvas[i]->Divide(2,1);
	Canvas[i]->cd(1);

	TString tmpString = TString(bnames[i]) + ">>H_" + TString(fnames[i]);
	std::cout << "tmpString = " << tmpString << std::endl; // "debug" printout
	hist[i] = new TH2F("H_"+TString(fnames[i]),TString(bnames[i]),100,range[k],range[k+1],100,range[k+2],range[k+3]);
	tmpTree->Draw(tmpString.Data(),"trackx@.GetEntries()==1 && tracky@.GetEntries()==1 && trackx.q>0 && tracky.q>0","colz");
	hist[i]->GetXaxis()->SetTitle(TString(title[m]));
	hist[i]->GetYaxis()->SetTitle(TString(title[m+1]));
	cmsprem->Draw();

	hprofile[i] = new TProfile("P_"+TString(fnames[i]),TString(bnames[i]),100,range[k],range[k+1],range[k+2],range[k+3]);
	tmpTree->Draw(TString(bnames[i])+">>P_"+TString(fnames[i]),"trackx@.GetEntries()==1 && tracky@.GetEntries()==1 && trackx.q>0 && tracky.q>0");
	hprofile[i]->GetXaxis()->SetTitle(TString(title[m]));
	hprofile[i]->GetYaxis()->SetTitle(TString(title[m+1]));
	cmsprem->Draw();

	for (int m=0; m<100; m++)
	for (int n=0; n<100; n++)
	{
		if (hist[i]->GetBinContent(m,n)>0) 
		{
			vecm.push_back(m); vecn.push_back(n);
		}
	}	
	fit_Low = hist[i]->GetXaxis()->GetBinCenter(vecm.front()-1);
	fit_High = hist[i]->GetXaxis()->GetBinCenter(vecm.back()+1);
	vecn.clear();	vecm.clear();		// clear the vector

	//fit_Low = 10.0;
	//fit_High = 80.0;
	hprofile[i]->GetXaxis()->SetRangeUser(fit_Low,fit_High);
//	hprofile[i]->SetMinimum(fit_Low);
	hprofile[i]->SetStats(0);
	hist[i]->SetStats(0);
	hist[i]->Draw();
	cmsprem->Draw();

	Canvas[i]->cd(2);
	hprofile[i]->Draw();
	cmsprem = new TLatex(-10,110,"CMS Preliminary");
	cmsprem->Draw();


	function[i] = new TF1("f1","[0]+([1]*x)",fit_Low,fit_High);
	//if (i==8) function[i] = new TF1("f1","expo",fit_Low,fit_High);

	function[i]->SetParameters(1,1);

	hprofile[i]->Fit("f1");
	hprofile[i]->Draw();
	cmsprem = new TLatex(0,108,"CMS Preliminary");
	cmsprem->Draw();

	TLegend *legend2 = new TLegend (0.14,0.95,.6,0.8,"Linear Fit");
	TF1 *fun=hprofile[i]->GetFunction("f1");
	cout<<"##############################################################"<<endl;
	cout<<"chi-sqr = "<<fun->GetChisquare()<<endl;
	cout<<"NDF = "<<fun->GetNDF()<<endl;
	o_file<<TString(fnames[i])<<"\t\t"<<setprecision(3)<<fixed<<fun->GetChisquare()/fun->GetNDF()<<"\t\t"<<fun->GetParameter(0)<<"\t\t"<<fun->GetParError(0)<<"\t\t\t"<<fun->GetParameter(1)<<"\t\t"<<fun->GetParError(1)<<endl; 
	cout<<"##############################################################"<<endl;
	sprintf(message,"Total : #chi^{2}/NDF = %.2d",fun->GetChisquare()/fun->GetNDF());
	legend2->AddEntry(fun,message);
	sprintf(message,"p0 = %.3f #pm  %.3f",fun->GetParameter(0),fun->GetParError(0));
	legend2->AddEntry(fun,message);
	sprintf(message,"p1 = %.3f  #pm %.3f",fun->GetParameter(1),fun->GetParError(1));
	legend2->AddEntry(fun,message);
	legend2->Draw();




//	Canvas[i]->SaveAs("Correlation_"+TString(fnames[i])+".png");
//	Canvas[i]->SaveAs("Correlation_"+TString(fnames[i])+".pdf");
	Canvas[i]->SaveAs(Form("Correlation_%s_Run%d.pdf",fnames[i],name));

}

o_file.close();
	}
	}