#include <iostream>
#include <cctype>
#include <cctype>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <utility>
#include "TChain.h"
#include "TH1.h"
#include "TPad.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TSystemDirectory.h"
#include "TList.h"
#include "TFile.h"
#include "/home/veronneau/RootUtils/AtlasUtils.h"

int btagging() {
	TChain *mchain = new TChain();
	double xbins[11]={20,30,40,50,60,75,90,110,140,200,14000};
	TH1D *trubint[4][8][3],*trub[4][8][3];
	TH1D *reco[3][4][8][3],*recoint[3][4][8][3];
	for (int m=0;m<3;++m) {
		for (int j=0;j<4;++j) {
			for (unsigned int l=0;l<8;++l) {
				trubint[j][l][m] = new TH1D(Form("trubint%d%d%d",j,l,m),"",10,xbins);
				trubint[j][l][m]->Sumw2();
				for (int i=0;i<3;++i) {
					recoint[i][j][l][m] = new TH1D(Form("recoint%d%d%d%d",i,j,l,m),"",10,xbins);
					recoint[i][j][l][m]->Sumw2();
				}
			}
		}	
	}
	TH1D *jetpt[8],*scaledjetpt[8],*scaledoverjetpt[8];
	for (int j=0;j<8;++j) {
		jetpt[j] = new TH1D(Form("jetpt%d",j),"jetpt",50,0.,500.);
		scaledjetpt[j] = new TH1D(Form("scaledjetpt%d",j),"scaledjetpt",50,0.,500.);
		scaledoverjetpt[j] = new TH1D(Form("scaledoverjetpt%d",j),"scaledoverjetpt",50,0.,500.);
		jetpt[j]->Sumw2();
		scaledjetpt[j]->Sumw2();
		scaledoverjetpt[j]->Sumw2();
	}
	TH1D *ratiojetpt[2];
	ratiojetpt[0] = new TH1D("ratiojetpt","ratiojetpt",50,0.,500.);
	ratiojetpt[1] = new TH1D("ratiooverjetpt","ratiooverjetpt",50,0.,500.);
	for (int i=0;i<2;++i) {ratiojetpt[i]->Sumw2();}
	TH1D *eff[3][4][8][3],*eff_mcatnlo[3][4][3],*sf[3][4][8][3];
	for (int m=0;m<3;++m) {
		for (int i=0;i<3;++i) {
			for (int j=0;j<4;++j) {
				eff_mcatnlo[i][j][m] = new TH1D(Form("eff_mcatnlo%d%d%d",i,j,m),"",10,xbins);
				eff_mcatnlo[i][j][m]->Sumw2();
				for (unsigned int l=0;l<8;++l) {
					eff[i][j][l][m] = new TH1D(Form("eff%d%d%d%d",i,j,l,m),"",10,xbins);
					sf[i][j][l][m] = new TH1D(Form("sf%d%d%d%d",i,j,l,m),"",10,xbins);
					eff[i][j][l][m]->Sumw2();
					sf[i][j][l][m]->Sumw2();
				}
			}	
		}	
	}
	string ds;
	string srcdir="/data/atlas/atlasdata3/arobic/SUSY/DGWH";
	string plotdir="/data/atlas/atlasdata3/arobic/SUSY/DGWH/macros/plots";
	string dirname,suffix;
	string outname="/data/atlas/atlasdata3/arobic/SUSY/DGWH/macros/histos/out";
	string dsname,substr,subsubstr,sample,filter;
	vector<string> inputfiles,dsescomp,dses,dsids,samples,samplesonly;
	unsigned int found;
	int op,nEntries;
	int nTotal=0;
	float fjetpt;
	int label;
	TList *files;
	TLegend *leg;
	Double_t xl1=.55, yl1=0.75, xl2=xl1+.4, yl2=yl1+.2;
	vector<string> sop;
	sop.push_back("60");
	sop.push_back("70");
	sop.push_back("80");
	vector<string> part;
	part.push_back("b");
	part.push_back("c");
	part.push_back("l");
	part.push_back("t");
	vector<int> scol;
	scol.push_back(kGreen);
	scol.push_back(kBlue);
	scol.push_back(kRed);
	//Branches
  //vector<float> *b_xsec;
  vector<double> *b_mcevtwgt;
  //vector<float> *b_sigjetpt;
  vector<float> *b_truthWpt;
  vector<float> *b_truthZpt;
	vector<int> *b_label;
	vector<double> *b_jeteta;
	vector<double> *b_jetphi;
	vector<double> *b_jetpt;
	vector<bool> *b_mv160;
	vector<bool> *b_mv170;
	vector<bool> *b_mv180;
	vector<int> *b_labelafterOR;
	vector<double> *b_jetetaafterOR;
	vector<double> *b_jetphiafterOR;
	vector<double> *b_jetptafterOR;
	vector<bool> *b_mv160afterOR;
	vector<bool> *b_mv170afterOR;
	vector<bool> *b_mv180afterOR;
	float xsec,lumibase,scale,nevents;
	unsigned pos;
	TFile *out;
	int cat=0;

	TCanvas *c1= new TCanvas("c1","c1",1500,500);
	TCanvas *c2= new TCanvas();
	TPad *pad1 = new TPad("pad1","pad1",0,0.27,1,1);
	TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.27);

	bool inclusive,last,nonWZ,cvetobveto,bfilter,massive;
	bool first;
	bool eta12, eta25;

	//Inclusive flag - to be set explicitely
	inclusive=false;
	massive=true;

	dirname="/data/atlas/atlasdata3/arobic/SUSY/DGWH/macros/histos";

	//Eta flags
	eta12=true;
	eta25=false;

	lumibase=20.3; //in fb - data lumi

	//Fill nevents map
	map<string,float> wgtmap;
	map<string,float>::iterator map_itr;
	float nevts;
	ifstream inFile;
	inFile.open(Form("%s/MC_evtwgt.txt",srcdir.c_str()));
	if (!inFile.is_open()) {
	  cerr << "Can't open input file " << Form("%s/MC_evtwgt.txt",srcdir.c_str()) << endl;
	  exit(1);
	}
	while (!inFile.eof()) {
  	inFile >> dsname;
  	inFile >> nevts;
		if (inFile.eof()) break; //Needed to avoid repeat of the last line
		found = dsname.find(string("#"));
  	if (found!=std::string::npos && found==0) continue;
		wgtmap.insert(std::make_pair(dsname,nevts));
	}
	inFile.close();

	//Read in MC samples
	inFile.open(Form("%s/MC_Alpgen.txt",srcdir.c_str()));
	if (!inFile.is_open()) {
	  cerr << "Can't open input file " << Form("%s/MC_Alpgen.txt",srcdir.c_str()) << endl;
	  exit(1);
	}
	while (!inFile.eof()) {
  	inFile >> dsname;
		if (inFile.eof()) break; //Needed to avoid repeat of the last line
		found = dsname.find(string("#"));
  	if (found!=std::string::npos && found==0) continue;
		inputfiles.push_back(dsname);
	}
	inFile.close();

	//Form ds string
	for (unsigned int f=0; f<inputfiles.size(); ++f) {
		stringstream stream(inputfiles.at(f));
		dsescomp.clear();
		int n=0;
		dsname="";
		while( getline(stream, substr, '.') )
		{
			if (n==3) {dsname+=substr; dsids.push_back(substr);}
			else if (n==6) dsname+=string("_")+substr;
			if (n==4) {
				int m=0;
				stringstream substream(substr);
				while( getline(substream, subsubstr, '_') )
				{
					if (strcmp(dsids.at(f).substr(0,2).c_str(),"15")==0) { 
						if (m==3) sample=subsubstr;
						else if ((m==4 && isdigit(subsubstr.substr(0,1).at(0)))) sample+=string("_")+subsubstr;
						else if ((m==4 && isalpha(subsubstr.substr(0,1).at(0)))) filter=subsubstr;
						else if ((m==5 && isalpha(subsubstr.substr(0,1).at(0)))) filter=subsubstr;
					} else {
						if (m==2) sample=subsubstr;
						else if ((m==3 && isdigit(subsubstr.substr(0,1).at(0)))) sample+=string("_")+subsubstr;
						else if ((m==3 && isalpha(subsubstr.substr(0,1).at(0)))) filter=subsubstr;
						else if ((m==4 && isalpha(subsubstr.substr(0,1).at(0)))) filter=subsubstr;
					}
					++m;
				}
			}
			++n;
		}
		if (!filter.empty()) filter=string("_")+filter;
		samplesonly.push_back(sample);
		samples.push_back(sample+filter);
		dses.push_back(dsname);
		sample.clear();
		filter.clear();
		dsname.clear();
	}

	TH1::AddDirectory(kFALSE);
 
	mchain->SetName("susytree");
	for (unsigned int f=0; f<dses.size(); ++f) {
		ds=dses.at(f);
		bfilter=false;
		first=false;
		if (samples.at(f).find("Np0") != string::npos) {bfilter=true; first=true;}
		cvetobveto=false;
		last=false;
		if (samples.at(f).find("Np3") != string::npos) {cvetobveto=true; last=true;}
		nonWZ=false;
		if (samples.at(f).find("W") == string::npos && samples.at(f).find("Z") == string::npos) nonWZ=true;
		cout << ds << " " << samples.at(f)<< endl;
		pos = samples.at(f).find("Np");
		if (pos != string::npos) {
			filter=samples.at(f).substr(1,pos-1);
			if (strcmp(filter.c_str(),"enu")==0) cat=1;
			else if (strcmp(filter.c_str(),"munu")==0) cat=2;
			else if (strcmp(filter.c_str(),"taunu")==0) cat=3;
			else if (strcmp(filter.c_str(),"bb")==0) cat=4;
		}

		//Filling chain by reading Root files	
		cout << dirname << endl;
		TSystemDirectory dir(dirname.c_str(), dirname.c_str());
		files = dir.GetListOfFiles();
		if (files) {
		   TSystemFile *file;
		   TString fname;
		   TIter next(files);
		   while ((file=(TSystemFile*)next())) {
		      fname = file->GetName();
		      if (!file->IsDirectory() && fname.EndsWith(".root") && fname.Contains(ds.c_str())) {
						mchain->Add(fname.Data());
		      }
		   }
		}
		
		//Set Branch Address
		mchain->SetBranchAddress("b_mcevtwgt", &b_mcevtwgt);
		mchain->SetBranchAddress("b_labelafterOR",&b_labelafterOR);
		mchain->SetBranchAddress("b_jetetaafterOR",&b_jetetaafterOR);
		mchain->SetBranchAddress("b_jetphiafterOR",&b_jetphiafterOR);
		mchain->SetBranchAddress("b_jetptafterOR",&b_jetptafterOR);
		mchain->SetBranchAddress("b_mv160afterOR",&b_mv160afterOR);
		mchain->SetBranchAddress("b_mv170afterOR",&b_mv170afterOR);
		mchain->SetBranchAddress("b_mv180afterOR",&b_mv180afterOR);
		//Save xsec
		nEntries = mchain->GetEntries();
		cout << "nEntries= " <<  nEntries << ", xsec= " << xsec << endl;
		//Clear histos
		for (int m=0;m<3;++m) {
			for (int j=0;j<4;++j) {
				for (unsigned int l=0;l<8;++l) {
					trubint[j][l][m]->Reset();
					for (int i=0;i<3;++i) {
						recoint[i][j][l][m]->Reset();
					}
				}
			}	
		}
		//Loop over tree to fill histos
		for (int j=0;j<nEntries; ++j) {
			mchain->GetEntry(j);
			cout << b_labelafterOR->size() << endl;
			if (j%10000==0) cout << j << " events done" << endl;
			for (unsigned int l=0;l<b_labelafterOR->size(); ++l) {
				if (fabs(b_jetetaafterOR->at(l)) < 2.5) {
					for (int m=0; m<4; ++m) {
						if (m==0) label=5; //b
						else if (m==1) label=4; //c
						else if (m==2) label=0; //light
						else if (m==3) label=15; //tau
						if (abs(b_labelafterOR->at(l))==label) {
							//eff vs eta
							trubint[m][cat][1]->Fill(b_jetetaafterOR->at(l),b_mcevtwgt->at(0));
							if (b_mv160afterOR->at(l)) recoint[0][m][cat][1]->Fill(b_jetetaafterOR->at(l),b_mcevtwgt->at(0));
							if (b_mv170afterOR->at(l)) recoint[1][m][cat][1]->Fill(b_jetetaafterOR->at(l),b_mcevtwgt->at(0));
							if (b_mv180afterOR->at(l)) recoint[2][m][cat][1]->Fill(b_jetetaafterOR->at(l),b_mcevtwgt->at(0));
							//eff vs pt
							if (fabs(b_jetetaafterOR->at(l)) < 1.2) {
								if (eta12) {
									//cout << b_jetptafterOR->at(l) <<endl;
									trubint[m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv160afterOR->at(l)) recoint[0][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv170afterOR->at(l)) recoint[1][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv180afterOR->at(l)) recoint[2][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
								}
							}
							else {
								if (eta25) {
									trubint[m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv160afterOR->at(l)) recoint[0][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv170afterOR->at(l)) recoint[1][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv180afterOR->at(l)) recoint[2][m][cat][0]->Fill(b_jetptafterOR->at(l),b_mcevtwgt->at(0));
								}
							}
							//eff vs phi
							if (fabs(b_jetetaafterOR->at(l)) < 1.2) {
								if (eta12) {
									//cout << b_jetptafterOR->at(l) <<endl;
									trubint[m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv160afterOR->at(l)) recoint[0][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv170afterOR->at(l)) recoint[1][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv180afterOR->at(l)) recoint[2][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
								}
							}
							else {
								if (eta25) {
									trubint[m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv160afterOR->at(l)) recoint[0][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv170afterOR->at(l)) recoint[1][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
									if (b_mv180afterOR->at(l)) recoint[2][m][cat][2]->Fill(b_jetphiafterOR->at(l),b_mcevtwgt->at(0));
								}
							}
						}
					}
				}
			}
		}

		mchain->Reset();

	}
	
}

void Book_TH1D(std::string HistoName, int nxbins, double xbins[], bool suffix)
{

	std::string suffixhisto;
  std::string eta12="_eta12";
  std::string eta25="_eta25";

	TH1D *hist;

	if (suffix) {
	  for (int i=0;i<2; ++i){
	    if (i==0) suffixhisto = HistoName+eta12;
	    else if (i==1) suffixhisto = HistoName+eta25;
	    else break;
	
			hist = new TH1D (suffixhisto.c_str(),suffixhisto.c_str(),nxbins,xbins);
 			hist->Sumw2();
 			th1dMap.insert(std::make_pair(suffixhisto,hist));
		}
	}
	else {
		hist = new TH1D (HistoName.c_str(),HistoName.c_str(),nxbins,xbins);
 		hist->Sumw2();
 		th1dMap.insert(std::make_pair(HistoName,hist));
	}

}

void Book_TH1D(std::string HistoName, int nxbins, float xmin, float xmax, bool suffix)
{

	std::string suffixhisto;
//  std::string sr="_sr";
//  std::string qcdvr="_qcdvr";
//  std::string ttbarcr="_ttbarcr";
//  std::string stwcr="_stwcr";
//  std::string vr1="_vr1";
//  std::string vr2="_vr2";
  std::string eta12="_eta12";
  std::string eta25="_eta25";

	TH1D *hist;

	if (suffix) {
	  for (int i=0;i<2; ++i){
//	    if (i==0) suffixhisto = HistoName+sr;
//	    else if (i==1) suffixhisto = HistoName+qcdvr;
//	    else if (i==2) suffixhisto = HistoName+ttbarcr;
//	    else if (i==3) suffixhisto = HistoName+stwcr;
//	    else if (i==4) suffixhisto = HistoName+vr1;
//	    else if (i==5) suffixhisto = HistoName+vr2;
	    if (i==0) suffixhisto = HistoName+eta12;
	    else if (i==1) suffixhisto = HistoName+eta25;
	    else break;
	
			hist = new TH1D (suffixhisto.c_str(),suffixhisto.c_str(),nxbins,xmin,xmax);
	  	hist->Sumw2();
	  	th1dMap.insert(std::make_pair(suffixhisto,hist));
		}
	}
	else {
		hist = new TH1D (HistoName.c_str(),HistoName.c_str(),nxbins,xmin,xmax);
  	hist->Sumw2();
  	th1dMap.insert(std::make_pair(HistoName,hist));
	}

}

void Fill_TH1D(std::string HistoName, float value_x, float weight, bool label, std::string binlabel)
{

  TH1D *hist=Retrieve_TH1D(HistoName);
  hist->Fill(value_x,weight);
  if (label) {
    int bin=(int)value_x+1;
    hist->GetXaxis()->SetBinLabel(bin,binlabel.c_str());
  }

}

void Fill_TH1D(std::string HistoName, float value_x, float weight, bool sr, bool qcdvr, bool stwcd, bool ttbarcr, bool vr1, bool vr2)
{

	std::string suffixhisto;
	
  for (int i=0;i<6; ++i){
		ok=true;
    if (i==0 && sr) suffixhisto = HistoName+string("_sr");
    else if (i==0 && !sr) ok=false;
    else if (i==1 && qcdvr) suffixhisto = HistoName+string("_qcdvr");
    else if (i==1 && !qcdvr) ok=false;
    else if (i==2 && ttbarcr) suffixhisto = HistoName+string("_ttbarcr");
    else if (i==2 && !ttbarcr) ok=false;
    else if (i==3 && stwcr) suffixhisto = HistoName+string("_stwcr");
    else if (i==3 && !stwcr) ok=false;
    else if (i==4 && vr1) suffixhisto = HistoName+string("_vr1");
    else if (i==4 && !vr1) ok=false;
    else if (i==5 && vr2) suffixhisto = HistoName+string("_vr2");
    else if (i==5 && !vr2) ok=false;
    else break;
		
		if (ok) {
			//cout << "suffixhisto= " << suffixhisto << endl;
  		TH1D *hist=Retrieve_TH1D(suffixhisto);
  		hist->Fill(value_x,weight);
		}
	}

}

TH1D* Retrieve_TH1D(std::string HistoName)
{

  std::map<std::string,TH1D *>::iterator itr = th1dMap.find(HistoName);

  // If the key can not be found
  if(itr==th1dMap.end()) {
    std::cerr << "Error: the histogram " << HistoName
          << " does not exist as a 1d histogram!" << std::endl;
    return 0;
  }

  return (*itr).second;

}