/*Reads in the lhco and weights.out files, before creating the log weight comparison -29/1/15 Ant
*/


#include <iostream>
#include <fstream>
#include <string>
#include <vector>

//root includes


#include "TROOT.h"
#include "TApplication.h"

#include "TPad.h"
#include "TChain.h"
#include "TString.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TFile.h"   
#include "THStack.h"
#include "TLegend.h"
#include "TStyle.h"
using namespace std;

//////////////////////////////////////////////////////
// Class to hold muon information read in from LHCO //
//////////////////////////////////////////////////////
class Muon{
	public:
		double phi,eta, Pt, m, charge;
		Muon (double PHI,double ETA,double PT,double M,double CHARGE):
			phi(PHI),eta(ETA),Pt(PT),m(M),charge(CHARGE){}
		Muon(){}
};

//////////////////////////////////////////////////////
// Class to hold jet information read in from LHCO //
//////////////////////////////////////////////////////
class Jet{
	public:
		double phi,eta,m, Pt;
		Jet(){}
		Jet(double PHI,double ETA,double PT,double M):
			phi(PHI),eta(ETA),m(M),Pt(PT){}
};

///////////////////////////////////////////////////////////
// Class to hold missingET information read in from LHCO //
///////////////////////////////////////////////////////////
class MissingET{
	public:
		double phi,eta,m;
		MissingET(){}
		MissingET(double PHI,double ETA,double M):
			phi(PHI),eta(ETA),m(M){}
};

///////////////////////////////////////////////////////////////
// Class to all information on an event in LHCO AND out file //
//
// Stores information on weight, error on weight, whether or 
// not the weight has been processed and all lines of the lhco
// file in the format we decided 
//
// Initalize using constructor - deletes it all 
///////////////////////////////////////////////////////////////
class Event{
	public:	
		int m_event_no;
		int m_run_no;
		Muon *m_mu1,*m_mu2;
		Jet  *m_j1,*m_j2,*m_j3,*m_j4;
		MissingET *m_mET;
		
		double m_weight,m_sig_weight;
		bool weighted;
	public:
		Event(){weighted = false;}
		Event(int event_no, int run_no,Muon* mu1,Muon* mu2,Jet *j1,Jet* j2, Jet* j3, Jet* j4,MissingET* missingET):
			m_event_no(event_no), m_run_no(run_no),
			m_mu1(mu1),m_mu2(mu2),
			m_j1(j1), m_j2(j2), m_j3(j3), m_j4(j4),
			m_mET(missingET){
				weighted = false; 
			}

		void SetWeight(double weight, double error){
			m_weight = weight; 
			m_sig_weight = error;
			weighted  = true;
		}
	
		//free up memory of all these damn pointers
		void Clean(){
			delete m_mu1; 
			delete m_mu2;
			delete m_j1; delete m_j2; delete m_j3; delete m_j4;
			delete m_mET;
		}
		
		//destructor - called when an object of type Event is deleted/ogoes out of scope - will auto free the memory
		~Event(){
			Clean();
		}
};
///////////////////////////////////////////////////////////////////////
//	A holder class for ALL events from LHCO and MEM file 	     //
//
// Stores a list of objects of type Event called Entries.
//
// Has two key methods : ReadLHCO and ReadWeights that set all data in 
// Entries to match those in file. Can then access the data via
// the square bracket operator, for example to access first muons
// Pt of the 1st read in evnet use:  [0]->m_mu1->Pt 
//
// NOTE: Use constructor which takes two paths, lhco and weights
// this will auto reads both files - will take a while!!
///////////////////////////////////////////////////////////////////////
class EventList{
	public:
	//VARIABLES
	//---------

		//A vector of Event class defined above to store all information about events
		vector<Event*> Entries;	//store the events
		
	//FUNCTIONS
	//---------

		//Provide a way to add another event to the list
		void Add(Event* ent){
			Entries.push_back(ent);
		}

		//Set the weigth(and error) of an event 
		void SetWeight(int event_no, double weight, double weight_error){
			bool set = false;			
			for(unsigned int i = 0; i < Entries.size(); i++){
				if(Entries[i]->m_event_no == event_no){
					Entries[i]->SetWeight(weight,weight_error);
					set = true;
				}
			}
		//	if(!set) cout << " failed to find event with event number " << event_no << endl;
			
		}	
		
		///////////////////////////////////////////////////////////
		// Read in LHCO file and fill Entries with event data 	 //
		///////////////////////////////////////////////////////////
		void ReadLHCO(const char* lhco_filepath){
			//Open the file
			//--------------------------------
			ifstream lhco_file(lhco_filepath);

			//Check to see if we could open these paths
			if(!lhco_file){
				cout << " FAILED TO OPEN LHCO FILE : " << lhco_filepath << endl;
				exit(1);
			}

			int i = 0;
			//Cycle over all information in the file until we reach the end of it
			while(lhco_file){
					if(i > 500) break; //over 500 = seg fault
					//if(i%100 == 0)cout << "processed lhco entry " << i << endl;
					
					i++;
					
					cout << "Processed lhco entry " << i << endl;
					//First Line - ie event number
					int line = -1, event_no = -1, run_no = -1;	//will reuse line variable as it is not going to be stored permanently! 
					lhco_file >> line;
					lhco_file >> event_no; 
					lhco_file >> run_no;

					//Muon 1
					int ID, p7,p8,p9,p10;
					double  eta, phi, Pt, m,charge;
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Muon* mu1=new Muon(eta,phi,Pt,m,charge);

					//muon 2
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Muon *mu2 = new Muon(eta,phi,Pt,m,charge);
					
					//jets
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Jet* j1 = new Jet(eta,phi,Pt,m);
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Jet* j2 = new Jet(eta,phi,Pt,m);
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Jet* j3 = new Jet(eta,phi,Pt,m);
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;
					Jet* j4 = new Jet(eta,phi,Pt,m);
					
					//,issing ET
					lhco_file >> line >> ID >> eta >> phi >> Pt >> m >> charge >> p7 >> p8 >> p9 >> p10;		
					MissingET* missing_ET =  new MissingET(phi,Pt,m);

					Add(new Event(event_no,run_no,mu1,mu2,j1,j2,j3,j4,missing_ET));
			}

			//finish reading the data so close the file
			//----------------------
			lhco_file.close();
		}
		///////////////////////////////////////////////////////////////////
		// Read in weights - note this should be called AFTER ReadLHCO	 //
		///////////////////////////////////////////////////////////////////
		void ReadWeights(const char* weights_filepath){
			//open the file  
			//-------------------------------------
			ifstream weights_file(weights_filepath);

			//check it actually opened
			if(!weights_file){
				cout << " FAILED TO OPEN WEIGHTS FILE : " << weights_filepath << endl;
				exit(1);
			}

			//cycle over all the weights until we reach end of file
			int i = 0;
			while(weights_file){
				if(i%100 == 0) cout << "processed weight entry " << i << endl;				

				int event_no, param_card, tf_card;
				double weight, weight_error;	
					
				weights_file >> event_no >> param_card >> tf_card >> weight >> weight_error;
				SetWeight(event_no,weight,weight_error);
			
				i++;
			}

			//finished reading things so close the file
			//----------------------------------------
			weights_file.close();
			

		}

		//Output entries that have weights 
		//-------------------------------------------------------------
		void Print(){
			for(unsigned int i = 0; i < Entries.size(); i++){
				if(Entries[i]->weighted)
					cout << Entries[i]->m_event_no << "\t" << Entries[i]->m_weight << endl;
			}
		}

		//Clean frees all memory held by the EntryList	- DON'T CALL unless you've read files in
		//---------------------------------------------------------------------------------------
		void Clean(){
			for(unsigned int i = 0; i < Entries.size();i++){
				delete Entries[i];
			}
		}

		//GetSize() returns number of entries that EventList currently has stored - both weighted and not!!
		//-------------------------------------------------------------------------------------------------
		unsigned int GetSize(){	return Entries.size();}
		unsigned int GetWeightedSize(){ 
			unsigned int size = 0;
			for(unsigned int i = 0 ; i <  GetSize();i++){
				if(Entries[i]->weighted)
					size++;
			}
			return size;
		}
		
		//Overwrite the [] operators to allow for event_list_object[0] return the 0th entry in the Entries vector
		//--------------------------------------------------------------------------------------------------------
		Event* operator[](unsigned int i){ return  Entries[i];}

		//Default constructor - left blank
		//--------------------------------	
		EventList(){}
		
		//////////////////////////////////////////////////////////////////////////////////////////
		// Constructor to open LHCO and Weights file, read in the data and be ready to be used  //
		//////////////////////////////////////////////////////////////////////////////////////////
		EventList(const char* lhco_filepath,const char *weights_filepath){
			cout << "Reading lhco file " << lhco_filepath << "...";
			ReadLHCO(lhco_filepath);
			cout << "found " << GetSize() << "entries."  << endl;
	
			cout << "Reading weight file " << weights_filepath << "...";
			ReadWeights(weights_filepath);
			cout << "found " << GetWeightedSize()  <<" entries. "<< endl;
			
			//cout << "Reading lhco file " << lhco_filepath << "...";
			//ReadLHCO(lhco_filepath);
			//cout << "found " << GetSize() << "entries." << endl;
		}
};

/////////////////////////////////////////////////////////////////////////////
//    Main Entry point of the progamme - what's called in the .sh file!    //
/////////////////////////////////////////////////////////////////////////////
void Process(){
	//declare a new object of type EventList to store all the lhco info & weights 	
	//read in the files given by the the .sh file 
	//TOCOMPARE:
	//
	//MEM ran with sig_1->2.0*sig_1 vs MEM ran with sig_1
	//Between 8TeV vs 13Tev -> TTH
	//TTH vs TTBB
	
	cout << "Starting .. "<< endl;
	EventList events_default("ttHfabrice_nobranches_030415.lhco","ttH_defaultTF_500weights.out" );
	EventList events_jet("ttHfabrice_nobranches_030415.lhco","ttH_jetTF_500weights.out");
	EventList events_jetmuon("ttHfabrice_nobranches_030415.lhco","ttH_jetmuonTF_500weights.out");
		


//cout << "starting second file processing" <<endl;
	//EventList events_2("ttHfabrice_nobranches_030415.lhco","ttH_defaultTF_040415_weights.out");
	


	///////////////////////////////////////////////////////////
	// 		Process the matched data		 //
	///////////////////////////////////////////////////////////

	//histograms to show stuff
	TH1F *histogram_default =new TH1F("Default MadWeight Transfer Function",";-log_{10}(weight_{defaultTF}) ;Number of events",25,25,50);
	TH1F *histogram_jet = new TH1F("Jet Transfer Function",";-log_{10}(weight_{jetTF}) ;Number of events",25,25,50);
	TH1F *histogram_jetmuon = new TH1F("Jet Muon Transfer Function",";-log_{10}(weight_{jetmuonTF}) ;Number of events",25,25,50);
	TH1F *histogram_jet_ratio = new TH1F("Jet TF / Default TF",";#frac{-log_{10}(weight_{jetTF})}{-log_{10}(weight_{defaultTF})} ;Number of events",25,0.9,1.1);
	TH1F *histogram_jetmuon_ratio = new TH1F("Jet Muon TF / Default TF",";#frac{-log_{10}(weight_{jetmuonTF})}{-log_{10}(weight_{defaultTF})} ;Number of events",25,0.9,1.1);
	

	histogram_default->Sumw2();
	histogram_jet->Sumw2();
	histogram_jetmuon->Sumw2();
	histogram_jet_ratio->Sumw2();
	histogram_jetmuon_ratio->Sumw2();

	
	//fill histogram(s)
	//--------------------------------------------------------
	for(unsigned int i = 0; i < events_default.GetSize();i++){
			cout << log10(events_default[i]->m_weight)<<endl;//log10(events_2[i]->m_weight) << endl;
			
			if(events_default[i]->weighted)
			if(events_default[i]->m_weight != 0.0)
			{
			histogram_default->Fill(-log10(events_default[i]->m_weight));///log10(events_2[i]->m_weight));
			}
		}

	for(unsigned int i = 0; i < events_jet.GetSize();i++){
			cout << log10(events_jet[i]->m_weight)<<endl;//log10(events_2[i]->m_weight) << endl;
						
			if(events_jet[i]->weighted)
			if(events_jet[i]->m_weight != 0.0)
			{
			histogram_jet->Fill(-log10(events_jet[i]->m_weight));///log10(events_2[i]->m_weight));
			}
		}
		
	for(unsigned int i = 0; i < events_jetmuon.GetSize();i++){
			cout << log10(events_jetmuon[i]->m_weight)<<endl;//log10(events_2[i]->m_weight) << endl;
			if(events_jetmuon[i]->weighted)
			if(events_jetmuon[i]->m_weight != 0.0)
			{			
			histogram_jetmuon->Fill(-log10(events_jetmuon[i]->m_weight));///log10(events_2[i]->m_weight));
			}		
		}
		
	for(unsigned int i = 0; i < events_jet.GetSize();i++){
			cout << log10(events_jet[i]->m_weight)/log10(events_default[i]->m_weight) << endl;
			if(events_jet[i]->weighted)
			if(events_jet[i]->m_weight != 0.0)
			if(events_default[i]->weighted)
			if(events_default[i]->m_weight != 0.0)
			{
			histogram_jet_ratio->Fill(-log10(events_jet[i]->m_weight)/-log10(events_default[i]->m_weight));
			}
		}
	for(unsigned int i = 0; i < events_jetmuon.GetSize();i++){
			cout << log10(events_jetmuon[i]->m_weight)/log10(events_default[i]->m_weight) << endl;
			if(events_jetmuon[i]->weighted)
			if(events_jetmuon[i]->m_weight != 0.0)
			if(events_default[i]->weighted)
			if(events_default[i]->m_weight != 0.0)
			{
			histogram_jetmuon_ratio->Fill(-log10(events_jetmuon[i]->m_weight)/-log10(events_default[i]->m_weight));
			}
		}
	///////////////////////////////////////////////////////////
	// 	Save the histogrmas to a root file 		 //
	///////////////////////////////////////////////////////////
	TFile f("blahblah.root", "recreate");
	histogram_default->Write();
	histogram_jet->Write();
	histogram_jetmuon->Write();
	histogram_jet_ratio->Write();
	histogram_jetmuon_ratio->Write();
	f.Close(); 


	///////////////////////////////////////////////////////////
	// 		Output processed histogram(s)		 //
	///////////////////////////////////////////////////////////
	gStyle->SetOptStat(1);
        TCanvas a1(" ","unnormalized weight",700,700);
	
	TPad *grid = new TPad("grid","",0,0,1,1);
	grid->Draw();
	grid->cd();
	grid->SetBottomMargin(0.3);

	


	grid->SetGrid();

	histogram_jet_ratio->Draw(""); //The first file printed always has too much white space for some reason. So this dummy file is printed to 						avoid the problem (The white space is caused by the grid definition above, so could just move that to later)
	a1.Clear();

	histogram_jet->Draw("E3 hist"); //the hist is needed to combine sumw2 and an error bar specification like E3
	histogram_jet->Draw("same");
	histogram_jet_ratio->GetYaxis()->SetLabelSize(0.02);
	histogram_jet_ratio->GetYaxis()->SetTitleSize(0.03);
	histogram_jet_ratio->GetXaxis()->SetLabelSize(0.02); //sets size of numbers on axis
	histogram_jet_ratio->GetXaxis()->SetTitleSize(0.02); //sets size of axis title
	histogram_jet_ratio->GetXaxis()->SetTitleOffset(2.0); //moves axis title
	//histogram_jet_ratio->Sumw2();
	a1.Print("ttH_jetvsdefault_weights_280415.eps");
	a1.Clear();
	
	histogram_jet->Draw("E3 hist"); //the hist is needed to combine sumw2 and an error bar specification like E3
	histogram_jet->Draw("same");
	histogram_jetmuon_ratio->GetYaxis()->SetLabelSize(0.02);
	histogram_jetmuon_ratio->GetYaxis()->SetTitleSize(0.03);
	histogram_jetmuon_ratio->GetXaxis()->SetLabelSize(0.02); //sets size of numbers on axis
	histogram_jetmuon_ratio->GetXaxis()->SetTitleSize(0.02); //sets size of axis title
	histogram_jetmuon_ratio->GetXaxis()->SetTitleOffset(2.0); //moves axis title
	//histogram_jetmuon_ratio->Sumw2(); //plots error bars on graph
	a1.Print("ttH_jetmuonvsdefault_weights_280415.eps");
	a1.Clear();
	

	//histogram_jet->Sumw2();
	histogram_jet->Draw("E3 hist"); //the hist is needed to combine sumw2 and an error bar specification like E3
	histogram_jet->Draw("same");
	a1.Print("ttH_jetTF_weights_280415.eps");
	a1.Clear();

	histogram_jetmuon->Draw("E3 hist"); //the hist is needed to combine sumw2 and an error bar specification like E3
	histogram_jetmuon->Draw("same");
	a1.Print("ttH_jetmuonTF_weights_280415.eps");
	a1.Clear();


	histogram_default->Draw("E3 hist"); //the hist is needed to combine sumw2 and an error bar specification like E3
	histogram_default->Draw("same");
	a1.Print("ttH_defaultTF_weights_280415.eps");
	a1.Clear();

	histogram_default->SetLineColor(gROOT->GetColor(2)->GetNumber());
	histogram_jet->SetLineColor(gROOT->GetColor(8)->GetNumber());
	histogram_jetmuon->SetLineStyle(3);
	histogram_jet->SetLineWidth(3);
	histogram_jetmuon->SetLineWidth(3);
	histogram_default->SetLineWidth(3);
	histogram_jet->SetLineWidth(2);	
	histogram_jetmuon->SetLineWidth(2);
	histogram_default->Draw(""); //the hist is needed to combine sumw2 and an error bar specification like E3	
	//histogram_jet->Draw("E3 hist");				
	//histogram_jetmuon->Draw("E3 hist");	
	histogram_default->Draw("same");
	histogram_jet->Draw("same");	
	histogram_jetmuon->Draw("same");

	gPad->Update();
	TPaveStats *st3 = (TPaveStats*)histogram_default->FindObject("stats"); //find stats box
	TPaveStats *st4 = (TPaveStats*)histogram_jet->FindObject("stats");
	TPaveStats *st5 = (TPaveStats*)histogram_jetmuon->FindObject("stats");
	histogram_default->(TH1::kNoStats);
	st3->InsertText("blahblah");
	st3->SetY1NDC(0.4); //sets new value for x starting position
	st3->SetY2NDC(0.5); //sets new value for x ending position
	st4->SetY1NDC(0.5);
	st4->SetY2NDC(0.6);
	st5->SetY1NDC(0.6);
	st5->SetY2NDC(0.7);


	leg = new TLegend(0.6,0.7,1.0,0.9);
	leg->AddEntry(histogram_default, "Default Transfer Function","l");	
	leg->AddEntry(histogram_jet, "Jet Transfer Function","l");	
	leg->AddEntry(histogram_jetmuon, "Jet Muon Transfer Function","l");	
	leg->SetTextSize(0.03);
	leg->Draw();
	//histogram_jet->Sumw2();  //Once you've done Sumw2() once on a histogram, dont need to do it again
	//histogram_jetmuon->Sumw2();
	//histogram_default->Sumw2();
	a1.Print("ttH_TFcomparison_weights_280415.eps");


	///////////////////////////////////////////////////////////
	// 	Clean up resources: close files &  free memory	 //
	///////////////////////////////////////////////////////////

	//EventList needs to have memory freed	
	events_default.Clean();
	events_jet.Clean();
	events_jetmuon.Clean();

	
	//free histogram memory
	delete histogram_default;
	delete histogram_jet;
	delete histogram_jetmuon;
	delete histogram_jet_ratio;
	delete histogram_jetmuon_ratio;
}