// This file will serve to read in bins and use the TMath::Sort()
// function to order the buns from highest to lowest and then re-
// draw every histogram with a 68% confidence interval and 95% 
// confidence interval layered on top
// --------------------------------------------------------------
// Kevin Ingles
// kingles@fnal.gov
// BinSorter.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("HistogramsAndCanvas_EOW_");
TString _root(".root");
TString saveDir("~/mypictures/ThesisPictures/with_noMIP/Sorted/");
Double_t X0(50.),Y0(50.),Width(600.),Height(600.);
bool setLog;
const TString fileType = ".pdf";
void BinSorter(TString histType, TString date, int SaveCanvas=0)
{
  // Hisotgram vectors: RMS and totalEnergyLoss
  std::vector<TH1*> RMSHistograms;
  // Canvas vectors: RMS and totalEnergyLoss
  std::vector<TCanvas*> RMSCanvases;
  // array vectors
  //std::vector<double*> histogramBins;
  //std::vector<int*> numberOfBins, indexBins;
  // For Neyman construction plot
  std::vector<TF1*> NeymanPlotFunctions;
  TF1 *func;
  if (histType.Contains("h_WireEnergyRMS"))
  {
	func = new TF1("Neyman","(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_WireEnergy"))
  {
	func = new TF1("Neyman","(x<0.012)*11+(x>0.012)*-3",0.0,0.0024);\
	func->SetTitle("Energy per wire; #Delta [GeV]; Log(p [GeV])");
	setLog = false;
  } */
  else if (histType.Contains("h_TotalEnergy"))
  {
	  func = new TF1("Neyman","(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("Neyman","(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("Neyman","(x<-5.0)*11+(x>-5.0)*4",-8.0,2.0);
	  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("Neyman","(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("Neyman","(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("Neyman","(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("Neyman","(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);
  NeymanPlotFunctions.push_back(func);

  std::cout << "Vectors created!\n";
  // Fill vector with file names
  std::vector<TString> histogramFiles;
  /*histogramFiles.push_back(histFiles+"0.2"+_root);
  histogramFiles.push_back(histFiles+"0.3"+_root);
  histogramFiles.push_back(histFiles+"0.4"+_root);
  histogramFiles.push_back(histFiles+"0.5"+_root);
  histogramFiles.push_back(histFiles+"0.6"+_root);
  histogramFiles.push_back(histFiles+"0.7"+_root);
  histogramFiles.push_back(histFiles+"0.8"+_root);
  histogramFiles.push_back(histFiles+"0.9"+_root);
  histogramFiles.push_back(histFiles+"1.0"+_root);
  histogramFiles.push_back(histFiles+"10.0"+_root);
  histogramFiles.push_back(histFiles+"50.0"+_root);*/
  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::vector<TString> histTypes;
  //histTypes.push_back("h_WireEnergyHistograms");
  //histTypes.push_back("h_WireEnergy");
  //histTypes.push_bacK("h_TotalEnergy");
  //histTypes.push_back("h_EMShowerCount");
  //for (auto const& histType : histTypes)
  //{
  // Array to store the RMS integral values
  //TString histType = "h_WireEnergyRMS";
  double integralValues[14];
  int i(0);
  // Open root files, clone RMS histograms and push them
  std::cout << "Entering for file loop to fill vectors.\n";
  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.\n";
      //TH1D *tempHist = (TH1D*) myHistogramFile.Get("h_RMSperEvent");
      if (histType.Contains("EMShower")) 
      {
	TH1I *workingHist2;
	workingHist2 = (TH1I*) myHistogramFile.Get(histType);
	workingHist2->Draw();
	gPad->Update();
	char test;
	scanf("%c",&test);
	workingHist2->SetName("h_"+momentum);
	workingHist2->SetDirectory(0);
	RMSHistograms.push_back(workingHist2);
	myHistogramFile.Close();
	TCanvas *workingCanvas;
	RMSCanvases.push_back(workingCanvas);
      }
      else 
      {
	TH1D *workingHist;
	workingHist = (TH1D*) myHistogramFile.Get(histType);workingHist = (TH1D*) myHistogramFile.Get(histType);
	workingHist->Draw();
	gPad->Update();
	char test;
	scanf("%c",&test);
	workingHist->SetName("h_"+momentum);
	workingHist->SetDirectory(0);
	RMSHistograms.push_back(workingHist);
	myHistogramFile.Close();
	TCanvas *workingCanvas;
	RMSCanvases.push_back(workingCanvas);
      }
    }
    
  TFile *SortedBins;  
  if (SaveCanvas)
  {
    SortedBins = new TFile("/pnfs/dune/persistent/users/kingles/Muon_Histograms/TreeReader_Histogram/with_noMIP/"+date+"/SortedHistograms/"+histType+"_SortedBins.root","NEW");
    if (SortedBins->IsOpen()) std::cout << "Output file is open.\n";
  }
  // Sort RMS histograms
  int j(0);  // Keep track of coinciding indices
  for (auto hist : RMSHistograms)
    {
      std::cout << "Inside for hist : RMSHistograms loop.\n";
      TH1D *copy = (TH1D*) hist->Clone("copy");
      std::cout << "Copy hist to 'copy'.\n";
    
      // Make associated sorting arrays
      int myBins = hist->GetNbinsX();
      double myArray[myBins];
      int sortArray[myBins];
      // Fill array with bin values 
      Double_t firstBinContent = hist->GetBinContent(0);
	firstBinContent += hist->GetBinContent(1);
	Double_t finalBinContent = hist->GetBinContent(myBins+1);
	finalBinContent += hist->GetBinContent(myBins);
	hist->SetBinContent(1,firstBinContent);
	hist->SetBinContent(myBins,finalBinContent);
      double integralValue = hist->Integral();
      for (int i = 0; i < myBins; i++)
	{
	  myArray[i] = hist->GetBinContent(i+1);
	}
      std::cout << "Filled myArray with hist bin values.\n";
    
      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++)
	{
	  int kthBin = sortArray[k];
	  if (first)
	    {
	      lowBin = kthBin;
	      highBin = kthBin;
	      first = false;
	    }
	  hSum += myArray[sortArray[k]];
	  double ratio = hSum / integralValue; 
	  if (ratio > 0.683) copy->SetBinContent(sortArray[k]+1,0);
	  else 
	    { 
	      if (kthBin < lowBin) lowBin = kthBin;
	      if (kthBin > highBin) highBin = kthBin;
	    }
	  //cout << myArray[sortArray[k]] << " " << hSum << " " << ratio << endl;
	  double temp = copy->GetBinContent(sortArray[k]+1);
	  //cout << "Copy bin content: " << temp << endl;
	}

      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("Log of the RMS energy per wire; Log(#sigma_{#bar{#Delta}}[GeV]); Entries");
	  else if (histType.Contains("h_EMShowerCount")) hist->SetTitle("Number of showers per event; Shower Count; Frequency");
	  else if (histType.Contains("h_RMS_div_Average")) hist->SetTitle("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("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("Log of RMS of non-MIP wires; _{noMIP}Log(#sigma_{#bar{#Delta}});Entries");
	  else if (histType.Contains("h_TotalEnergy")) hist->SetTitle("Log of total energy deposited per event;Log(#Delta_{tot}[GeV]); Entries");
	  else if (histType.Contains("h_Average_noMIP")) hist->SetTitle("Log of average energy of non-MIP wires; _{noMIP}Log(#bar{#Delta}[GeV]);Entries");
	  else if (histType.Contains("h_WireEnergyAverage")) hist->SetTitle("Log of Average energy on wire; Log(#bar{#Delta}[GeV]);Entries");
    
      gStyle->SetOptStat("nemruoi");
      hist->Draw();
      copy->SetFillColor(17);
      copy->SetAxisRange(lowerLimit,upperLimit);
      copy->Draw("SAME");
      std::cout << "Write " << momentum << " canvas to TFile.\n";
      if (SaveCanvas == 1) 
	{
	  RMSCanvases[j]->Write();
	  RMSCanvases[j]->SaveAs(saveDir+histType+"_"+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();
      NeymanPlotFunctions.push_back(tempFunc);
    }

  TCanvas *NeymanPlot = new TCanvas("c_NeymanPlot","NeymanPlot",X0,Y0,Width,Height);
  bool firstDraw = true;
  for (auto plot : NeymanPlotFunctions)
    {
      if (firstDraw) 
	{
	  plot->Draw();
	  firstDraw = false; 
	}
      plot->Draw("same");
    }
  if (SaveCanvas == 1) 
    {
      NeymanPlot->Write();
      NeymanPlot->SaveAs(saveDir+"Neyman_"+histType+date+fileType);
    }
   if (SaveCanvas) SortedBins->Close();
  //}
}