#include <vector>
#include <iostream>
#include <fstream>
#include "TFile.h"
#include "TLine.h"
#include "TEllipse.h"
#include "TGraphErrors.h"
#include "TGAxis.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TPad.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TPolyLine.h"
#include "TProfile.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TLatex.h"
#include "TPaveStats.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TFile.h"
#include "TMath.h"
#include "THStack.h"
#include "TSystem.h"

#include "Util.cpp"

using namespace Utils;

void BunchTrainsOverlay(Int_t NOverlay,
                        TString outFileName,
                        vector<TString> fileNames, 
                        vector<TString> prefix, 
                        vector<TString> suffix, 
                        vector<TString> tags,
                        vector<TString> plotDir, 
                        vector<double> weights,
                        TString outDir = "./", 
                        TString psFile = "TowerDPD_JetRec_JetValidation",
                        TString plotType = ".pdf",
                        bool dataMCCompare = true,
                        bool dataMCRatio = true,
                        bool drawNorm    = true,
                        TString newLabel = "",
                        bool writeFile = true,
                        bool doLogY = false,
                        bool doLogZ = false) {

  cout << "----------------------------------------------" << endl;
  cout << "Starting plot overlay for bunch trains" << endl;
  cout << "   N overlay = " << NOverlay << endl;
  cout << "----------------------------------------------" << endl;
  
  ///-----------------------
  /// Output
  ///-----------------------
  TFile *outfile = new TFile(TString(outDir + outFileName),"RECREATE");
  
  ///-----------------------
  /// Input
  ///-----------------------
  vector<TFile*> files;   
  for (Int_t fileItr=0; fileItr < NOverlay; fileItr++) {
    cout << "Adding file: " << fileNames.at(fileItr) << endl;
    TFile *f = new TFile(fileNames.at(fileItr),"READ");
    if (f->IsOpen()) files.push_back(f);
    else cout << "Problem opening file: " << fileNames.at(fileItr) << endl;
  }
  cout << "Done Adding files. Size = " << files.size() << endl;
  cout << "Validity check: #1" << (files.at(0))->GetName() << endl;
 
  ///-----------------------
  /// Starting output PS
  ///-----------------------
  // Canvas
  TCanvas *canvas = (TCanvas*)gROOT->FindObject("c1");
  if (canvas) {canvas->Clear();} 
  else {
    style();
    canvas = (TCanvas*)gROOT->FindObject("c1");
  }
  canvas->Print(TString(psFile + ".ps["));
  TGaxis::SetMaxDigits(3);
  
  ///-----------------------
  /// Plot list
  ///-----------------------
  
  const int nBCIDbins = 9;
  const int nEtaBins  = 8;  
  Float_t etaBins[9]  = {0.0,0.3,0.8,1.2,1.7,2.1,2.8,3.6,4.5};
  
  Int_t n1DPlots = 14;
  TString plot1DList[50] = {
                            "tower_e",
                            "tower_et",
                            
                            "em_clus_pt",
                            "em_clus_pt",
                            "em_clus_isolation",
                            "em_clus_first_eng_dens",
                            "em_clus_lateral",
                            "em_clus_longitudinal",
                            "em_clus_center_lambda",
                            
                            "lc_clus_e",
                            "lc_clus_pt",
                            
                            "jet_e",
                            "jet_pt",
                            "jet_width"
                           };

  ///-----------------------                      
  /// 1D Plots   
  ///-----------------------
  cout << "Plotting 1D plots" << endl;
  for (Int_t i=0; i<n1DPlots; i++) {
    TH1F* plots[50];
    Float_t g_mean_pts[8];
    Float_t g_mean_error[8];
    Float_t g_RMS_pts[8];
    for (Int_t etaItr = 0; etaItr < nEtaBins; etaItr++) {
      for (Int_t bcidItr=0; bcidItr < nBCIDbins; bcidItr++) {
        cout << "Reading bcidItr = " << bcidItr << " and plot #" << i << endl;
        TH1F *h = (TH1F*)((files.at(bcidItr))->Get(TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i] + Form("_BCID%d_%d_to_%d",bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)))));
        if (h) cout << "Got histogram: " << TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i] + Form("_BCID%d_%d_to_%d",bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))) << endl;
        else { 
          cout << "Did not get histogram: " << TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i] + Form("_BCID%d_%d_to_%d",bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))) << endl;
          cout << "Trying to get vertex histogram without suffix" << endl;
          h = (TH1F*)((files.at(bcidItr))->Get(TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i])));
          if (h) cout << "Got vertex histogram: " << TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i]) << endl;
          else cout << "Still did not get histogram: " << TString(plotDir.at(bcidItr) + prefix.at(bcidItr) + plot1DList[i]) << endl;
        }
        plots[bcidItr] = h;
        g_mean_pts[bcidItr] = plots[bcidItr]->GetMean();
        g_mean_error[bcidItr] = plots[bcidItr]->GetMeanError();
        g_RMS_pts[bcidItr] = plots[bcidItr]->GetRMS();
      }
      
      Float_t YMin = 0.1; Float_t YMax = 0.0;
      bool logx = kFALSE;
      bool logy = kFALSE;
      bool doFit = kFALSE;
      bool drawStats = kFALSE;
      bool drawExtra = kTRUE;
      
      if (logy && (YMax < plots[0]->GetMaximum()*10000.)) YMax = plots[0]->GetMaximum()*10000.;
      else if (YMax < plots[0]->GetMaximum()*1.5) YMax = plots[0]->GetMaximum()*1.5;
  
      Float_t gBCID[8] = {1,2,3,4,5,6,7,8};
      Float_t gBCIDerror[8] = {0,0,0,0,0,0,0,0};
      TGraphErrors *g_mean = new TGraphErrors(8,gBCID,g_mean_pts,gBCIDerror,g_mean_error);
      TGraphErrors *g_RMS = new TGraphErrors(8,gBCID,g_mean_pts,gBCIDerror,g_RMS_pts);
      TCanvas *canv_graph = new TCanvas(TString(plot1DList[i] + Form("_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))),TString(plot1DList[i] + Form("_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))));
      g_mean->SetMarkerSize(1.4); g_mean->SetLineColor(kBlack);
      g_RMS->SetMarkerSize(0.2);  g_RMS->SetLineColor(kRed);    g_RMS->SetLineWidth(1.3);
      g_RMS->Draw("ap");
      g_mean->Draw("p");
      g_RMS->GetXaxis()->SetTitle("Position within the bunch train");
      g_RMS->GetYaxis()->SetTitle(plots[0]->GetXaxis()->GetTitle());
      outfile->cd();
      g_mean->Write(TString(plot1DList[i] + Form("_mean_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))));
      g_RMS->Write(TString(plot1DList[i] + Form("_RMS_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))));
      canv_graph->Write(TString(plot1DList[i] + Form("_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))),TObject::kOverwrite);
  
      canvas =
          PlotAll(plots,
                  nBCIDbins,
                  tags,
                  "",
                  "",
                  (const char*)(outDir + "overlay_" + plot1DList[i] + Form("_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)) + plotType),
                  (const char*)("overlay_" + plot1DList[i] + Form("_%d_to_%d",(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.))),
                  outfile,
                  0.0, 0.0,
                  YMin, YMax,
                  drawNorm,
                  logx, logy,
                  drawStats,
                  drawExtra,
                  writeFile,
                  dataMCCompare,
                  dataMCRatio,
                  newLabel
                 );
      canvas->Print(TString(psFile + ".ps"));
    }
  }
  cout << "Finished Plotting 1D plots" << endl;

/*
  ///------------------
  /// Set histograms
  ///------------------
  for (Int_t bcidItr = 0; bcidItr < 9; bcidItr++) {
    for (Int_t etaItr = 0; etaItr < nEtaBins; etaItr++) {
      
      // Towers
      h_tower_e[bcidItr][etaItr]  = new TH1F(Form("tower_e_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Tower energy [MeV]; Number of towers: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nTowerEBins,minTowerE,maxTowerE);
      h_tower_et[bcidItr][etaItr] = new TH1F(Form("tower_et_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Tower E_{T} [MeV]; Number of towers: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nTowerEBins,minTowerE,maxTowerE);

      // EM Clusters
      h_em_clus_e[bcidItr][etaItr]               = new TH1F(Form("em_clus_e_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster energy [GeV]; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nClusterEBins,minClusterE,maxClusterE);
      h_em_clus_pt[bcidItr][etaItr]              = new TH1F(Form("em_clus_pt_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster p_{T} [GeV]; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nClusterEBins,minClusterE,maxClusterE);
      h_em_clus_isolation[bcidItr][etaItr]       = new TH1F(Form("em_clus_isolation_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster isolation; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),100.,0.0,1.0);
      h_em_clus_first_eng_dens[bcidItr][etaItr]  = new TH1F(Form("em_clus_first_eng_dens_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Normalized EM cluster energy density; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),200.,0.0,0.01);
      h_em_clus_lateral[bcidItr][etaItr]         = new TH1F(Form("em_clus_lateral_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster lateral moment; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),200.,0.0,1.0);
      h_em_clus_longitudinal[bcidItr][etaItr]    = new TH1F(Form("em_clus_longitudinal_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster longitudinal moment; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),200.,0.0,1.0);
      h_em_clus_center_lambda[bcidItr][etaItr]   = new TH1F(Form("em_clus_center_lambda_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";EM Cluster center lambda; Number of EM clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),500.,0.,1000.);

      // LC clusters
      h_lc_clus_e[bcidItr][etaItr]  = new TH1F(Form("lc_clus_e_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";LC Cluster energy [GeV]; Number of LC clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nClusterEBins,minClusterE,maxClusterE);
      h_lc_clus_pt[bcidItr][etaItr] = new TH1F(Form("lc_clus_pt_BCID%d_%d_to_%d", bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";LC Cluster p_{T} [GeV]; Number of LC clusters: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nClusterEBins,minClusterE,maxClusterE);
      
      // Jets
      h_jet_e     [bcidItr][etaItr]  = new TH1F(Form("jet_e_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Jet energy [GeV]; Number of jets: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nJetEBins,minJetE,maxJetE);
      h_jet_pt    [bcidItr][etaItr]  = new TH1F(Form("jet_pt_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Jet E_{T} [GeV]; Number of jets: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),nJetEBins,minJetE,maxJetE);
      h_jet_width [bcidItr][etaItr]  = new TH1F(Form("jet_width_BCID%d_%d_to_%d",  bcidItr,(Int_t)(etaBins[etaItr]*10.),(Int_t)(etaBins[etaItr+1]*10.)),Form(";Jet width; Number of jets: %1.2f < |#eta| < %1.2f",etaBins[etaItr],etaBins[etaItr+1]),40,0,0.8);
    }
  }
*/

}