/********************************
 * Kevin Ingles
 * kingles@fnal.gov
 * BinSorter_FC.C
 ********************************
 */
// ROOT 
#include "TH1.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLatex.h"
// STL
#include <iostream>
#include <vector>
// prototypes
// global variables
TString histFiles("FCHistogramsAndCanvas_EOW_");
TString feldcous("FeldmanCousins");
TString _root(".root");
TString saveDir("~/mypictures/ThesisPictures/with_noMIP/Sorted/");
TString fileType(".pdf");
Double_t X0(50.),Y0(50.),Width(600.),Height(600.);
bool setLog;
void BinSorter(TString histType, TString date, int SaveCanvas=0)
{
	std::cout << "Inside programs\n";
	// Hisotgram vectors: RMS and totalEnergyLoss
	std::vector<TH1*> RMSHistograms;
	// Canvas vectors: RMS and totalEnergyLoss
	std::vector<TCanvas*> RMSCanvases;
	std::vector<TF1*> FCPlotFunctions;
	TF1 *func;
	if (histType.Contains("h_WireEnergyRMS"))
	{
		func = new TF1("FC","(x<-5.0)*11+(x>-5.0)*4",-8.0,2.0);
		func->SetTitle("Log of RMS of energy on a wire; Log(#sigma_{#Delta}[GeV]); Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_TotalEnergy"))
	{
		func = new TF1("FC","(x<TMath::Log(50))*11+(x>TMath::Log(50))*4",0,10);
		func->SetTitle("Log of total energy deposited; Log(#Delta_{tot} [GeV]); Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_EMShowerCount"))
	{
		func = new TF1("FC","(x<5.0)*11+(x>5.0)*4",0,30);
		func->SetTitle("EM Showers; Number of Showers; Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_WireEnergyAverage"))
	{
		//func = new TF1("FC","(x<-5.0)*11+(x>-5.0)*4",-8.0,0.0); ///< original
		//func = new TF1("FC","(x<0.0)*11 + (x>0.0)*4",TMath::ATan(-1.0),TMath::ATan(0.0)); ///< original transformed
		func = new TF1("FC","(x<-5.0)*11+(x>-5.0)*4",-8.0,2.0); ///< extended
		//func = new TF1("FC","(x<0.0)*11 + (x>0.0)*4",TMath::ATan(-1.0),TMath::ATan(1.0/8.0)); ///< extended transformed
		func->SetTitle("Average energy on a wire; #bar{#Delta} [GeV]; Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_RMS_div_Average"))
	{
		func = new TF1("FC","(x<0.5)*11+(x>0.5)*4",-1.0,2.0);
		func->SetTitle("Log of RMS divided by average for energy on a wire; Log(#sigma_{#Delta}/#bar{#Delta}); Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_Average_noMIP"))
	{
		func = new TF1("FC","(x<-5.0)*11+(x>-5.0)*4",-8.0,1.0);
		func->SetTitle("Average energy on non-MIP wires; #Delta [GeV]; Log(p [GeV])");
		setLog = false;
	}
	else if (histType.Contains("h_RMS_noMIP"))
	{
		func = new TF1("FC","(x<-4.0)*11+(x>-4.0)*4",-8.0,3.0);
		func->SetTitle("Log of RMS of non-MIP wires; Log(#sigma_{#Delta} [GeV]); Log(p [GeV])");
		setLog = false;
	}
    else if (histType.Contains("h_RMS_div_Average_noMIP"))
	{
		func = new TF1("FC","(x<0.5)*11+(x>0.5)*4",-1.0,2.0);
		func->SetTitle("Log of RMS divided by average for energy on non-MIP wires; Log(#sigma_{#Delta}/#bar{#Delta}); Log(p [GeV])");
		setLog = false;
	}
	func->SetLineColor(kWhite);
	FCPlotFunctions.push_back(func);
	
	
	// Fill vector with file names
	std::vector<TString> histogramFiles;
	histogramFiles.push_back(histFiles+"100.0"+_root);
	histogramFiles.push_back(histFiles+"500.0"+_root);
	histogramFiles.push_back(histFiles+"1000.0"+_root);
	histogramFiles.push_back(histFiles+"5000.0"+_root);
	histogramFiles.push_back(histFiles+"10000.0"+_root);
	histogramFiles.push_back(histFiles+"20000.0"+_root);
	histogramFiles.push_back(histFiles+"50000.0"+_root);
	
	std::cout << "Making histogram vector\n";
	int q(0);
	for (auto file : histogramFiles)
    {
		// Make cloned histogram
		TString momentum = file;
		momentum.ReplaceAll(histFiles,"EOW_");
		momentum.ReplaceAll(_root,"");
		TFile myHistogramFile(file,"r");
		if (myHistogramFile.IsOpen()) std::cout << "File is open.  " << q << endl;
		TH1 *workingHist;
		if (histType.Contains("EMShower")) workingHist = (TH1I*) myHistogramFile.Get(histType);
		else workingHist = (TH1D*) myHistogramFile.Get(histType);
		if (workingHist == NULL) printf("Histogram pointer is null pointer.\n");
		else printf("Histogram pointer is valid.\n");
		workingHist->Draw("HIST"); //gPad->WaitPrimitive();
		gPad->Update();
		char test;
		printf("Debugging pause. Please press a character then eneter: \n");
		scanf("%c",&test);
		workingHist->SetName("h_"+momentum);
		workingHist->SetDirectory(0);
		
		// Set first and last in contents and normalize
		Int_t bins = workingHist->GetNbinsX();
		Double_t first, second, penult, final;
		first  = workingHist->GetBinContent(0);
		second = workingHist->GetBinContent(1);
		penult = workingHist->GetBinContent(bins);
		final  = workingHist->GetBinContent(bins+1);
		workingHist->SetBinContent(1,first+second);
		workingHist->SetBinContent(bins,penult+final);
		workingHist->SetBinContent(0,0);
		workingHist->SetBinContent(bins+1,0);
		workingHist->Smooth();
		Double_t norm = 1.0/workingHist->Integral();
		printf("The normalization constant is: %lf\n",norm);
		workingHist->Scale(norm);
		
		RMSHistograms.push_back(workingHist);
		myHistogramFile.Close();
		TCanvas *workingCanvas;
		RMSCanvases.push_back(workingCanvas);
		q++;
    }
	
	TFile *SortedBins;
	if (SaveCanvas)
	{
	  SortedBins = new TFile("/pnfs/dune/persistent/users/kingles/Muon_Histograms/TreeReader_Histogram/with_noMIP/"+date+"/SortedHistograms/"+histType+"_"+feldcous+"_SortedBins.root","NEW");
	  if (!SortedBins->IsOpen()) { printf("Failed to open TFile.\n"); exit(-1); }
	}
	
	std::cout << "Looping over histograms\n";
	int j(0);
	for (auto hist : RMSHistograms)
    {
		TH1D *copy = (TH1D*) hist->Clone("copy");
		hist->Draw("HIST");
		gPad()->Update();
				
		// Make associated sorting arrays
		int myBins = hist->GetNbinsX();
		/*Double_t firstBinContent = hist->GetBinContent(0);
		Double_t secondBinContent = hist->GetBinContent(1);
		Double_t finalBinContent = hist->GetBinContent(myBins+1);
		Double_t penultimateBinContent = hist->GetBinContent(myBins);
		hist->SetBinContent(1,firstBinContent+secondBinContent);
		hist->SetBinContent(myBins,finalBinContent+penultimateBinContent);
		Double_t norm1 = 1.0/hist->Integral();
		hist->Scale(norm1);
		Double_t norm2 = 1.0/copy->Integral();
		copy->Scale(norm2);
		
		//hist->SetBinContent(1,secondBinContent*norm1);
		//hist->SetBinContent(myBins,penultimateBinContent*norm1);*/
		Double_t maxBin(0.0); //= hist->GetBinContent(hist->GetMaximumBin());
		double myArray[myBins];
		int sortArray[myBins];
		// Fill array with bin values
		for (int i = 0; i < myBins; i++)
		{
			maxBin = -1.0;
			for (auto const& tempHist : RMSHistograms)
			{
				double tempMax = tempHist->GetBinContent(i+1);
				if (tempMax > maxBin) maxBin = tempMax;
			}
			if (maxBin <= 0.0) 
				myArray[i] = hist->GetBinContent(i+1);
			else
				myArray[i] = hist->GetBinContent(i+1)/maxBin;
		}
		
		TMath::Sort((Int_t)myBins,myArray,sortArray);
		
		double hsum(0);
		bool first = true;
		int lowBin(0), highBin(0);
		for (int k = 0; k <= myBins; k++)
		{
			// We include the '+1' since we are using the 
			// the entries directly from the histogram and
			// not myArray. Which means bins start 1.
			int kthBin = sortArray[k];
			if (first)
			{
				lowBin = kthBin;
				highBin = kthBin;
				first = false;
			}
			hsum += hist->GetBinContent(sortArray[k]+1); 
			printf("%3d \t %3d \t %10.9e\n \t %10.9e\n", k, sortArray[k], hist->GetBinContent(sortArray[k]+1), hsum);
			//std::cout << hsum << endl;
			if (hsum > 0.683) copy->SetBinContent(sortArray[k]+1,0);
			else 
			{ 
				if (kthBin < lowBin) lowBin = kthBin;
				if (kthBin > highBin) highBin = kthBin;
			}
		}
		
		Double_t lowerLimit = hist->GetBinCenter(lowBin + 1);
		Double_t upperLimit = hist->GetBinCenter(highBin + 1);

		TString histName = (TString) hist->GetName();
		TString momentum = histName;
		momentum.ReplaceAll(histFiles,"");
		RMSCanvases[j] = new TCanvas("c_"+momentum,momentum,X0,Y0,Width,Height);
		X0 += 10; Y0 += 10; 
    
		if (histType.Contains("h_WireEnergyRMS")) hist->SetTitle("FC Log of the RMS energy per wire; Log(#sigma_{#bar{#Delta}}[GeV]); Entries");
		else if (histType.Contains("h_EMShowerCount")) hist->SetTitle("FC Number of showers per event; Shower Count; Frequency");
		else if (histType.Contains("h_RMS_div_Average")) hist->SetTitle("FC Log of RMS of energy on wire divided by average energy on wire;Log(#sigma_{#bar{#Delta}}/#bar{#Delta});Entries");
		else if (histType.Contains("h_RMS_div_Average_noMIP")) hist->SetTitle("FC Log of RMS of non-MIP wires divided by average energy on non-MIP wire; _{noMIP}Log(#sigma_{#bar{Delta}}/#bar{#Delta});Entries");
		else if (histType.Contains("h_RMS_noMIP")) hist->SetTitle("FC Log of RMS of non-MIP wires; _{noMIP}Log(#sigma_{#bar{#Delta}});Entries");
		else if (histType.Contains("h_TotalEnergy")) hist->SetTitle("FC Log of total energy deposited per event;Log(#Delta_{tot}[GeV]); Entries");
		else if (histType.Contains("h_Average_noMIP")) hist->SetTitle("FC Log of average energy of non-MIP wires; _{noMIP}Log(#bar{#Delta}[GeV]);Entries");
		else if (histType.Contains("h_WireEnergyAverage")) hist->SetTitle("FC Log of Average energy on wire; Log(#bar{#Delta}[GeV]);Entries");
	
		gStyle->SetOptStat("nemruoi");
		hist->Draw("HIST");
		copy->SetFillColor(17);
		copy->SetAxisRange(lowerLimit,upperLimit);
		copy->Draw("HIST SAME");
		std::cout << "Write " << momentum << " canvas to TFile.\n";
		if (SaveCanvas == 1) 
		{
			RMSCanvases[j]->Write();
			RMSCanvases[j]->SaveAs(saveDir+histType+"_"+feldcous+"_"+momentum+date+fileType);
		}
		j++;
    
		if (setLog)
		{
			lowerLimit = TMath::Log(lowerLimit);
			upperLimit = TMath::Log(upperLimit);
		}
		TString energy = momentum;
		energy.ReplaceAll("h_EOW_","");
		TF1 *tempFunc = new TF1(momentum,"TMath::Log("+energy+")",lowerLimit,upperLimit);
	  
		if (SaveCanvas == 1) tempFunc->Write();
		FCPlotFunctions.push_back(tempFunc);
    }

	std::cout << "Drawing final canvases\n";
	TCanvas *FCPlot = new TCanvas("c_FeldmanCousins","FeldmanCousins",X0,Y0,Width,Height);
	bool firstDraw = true;
	for (auto plot : FCPlotFunctions)
    {
		if (firstDraw) 
		{
			plot->Draw();
			firstDraw = false; 
		}
		else plot->Draw("same");
    }
	if (SaveCanvas == 1) 
    {
      FCPlot->Write();
      FCPlot->SaveAs(saveDir+feldcous+"_"+histType+date+fileType);
    }
  if (SaveCanvas) SortedBins->Close();
}