#include <iostream>
#include "TString.h"
#include "TFile.h"
#include "TH1F.h"
#include "TKey.h"
#include "TROOT.h"
#include "TClass.h"
#include "TMath.h"
#include "TGraph.h"
#include <vector>
#include <map>
using namespace std;

int n_SR = 7;
TString SR[] =  {"ATLAS_Bkg_LH_Other" ,  "ATLAS_Bkg_LH_TT" ,  "ATLAS_Bkg_LH_Wlnu" ,  
		 "ATLAS_Bkg_LH_Zleplep" ,  "ATLAS_Bkg_LH_Ztautau" ,  "ATLAS_Bkg_LH_qcd_cC"};

int n_CR = 6;
TString CR_D[] = {"ATLAS_Bkg_LH_Other_cD", "ATLAS_Bkg_LH_TT_cD", "ATLAS_Bkg_LH_Wlnu_cD", 
		  "ATLAS_Bkg_LH_Zleplep_cD", "ATLAS_Bkg_LH_Ztautau_cD"};
TString CR_B[] = {"ATLAS_Bkg_LH_Other_cB", "ATLAS_Bkg_LH_TT_cB", "ATLAS_Bkg_LH_Wlnu_cB", 
		  "ATLAS_Bkg_LH_Zleplep_cB", "ATLAS_Bkg_LH_Ztautau_cB"};

double other_sr[100];
double tt_sr[100];
double wlnu_sr[100];
double zleplep_sr[100];
double ztautau_sr[100];
double qcd_sr[100];

double other_sr_err[100];
double tt_sr_err[100];
double wlnu_sr_err[100];
double zleplep_sr_err[100];
double ztautau_sr_err[100];
double qcd_sr_err[100];


double other_b[100];
double tt_b[100];
double wlnu_b[100];
double zleplep_b[100];
double ztautau_b[100];
double qcd_b[100];

double other_d[100];
double tt_d[100];
double wlnu_d[100];
double zleplep_d[100];
double ztautau_d[100];
double qcd_d[100];


// add some signal samples for tests
double mA100_b[100], mA100_g[100], mA100_be[100], mA100_ge[100];
double mA200_b[100], mA200_g[100], mA200_be[100], mA200_ge[100];
double mA300_b[100], mA300_g[100], mA300_be[100], mA300_ge[100];
double mA400_b[100], mA400_g[100], mA400_be[100], mA400_ge[100];
double mA500_b[100], mA500_g[100], mA500_be[100], mA500_ge[100];
double mA600_b[100], mA600_g[100], mA600_be[100], mA600_ge[100];
double mA700_b[100], mA700_g[100], mA700_be[100], mA700_ge[100];
double mA800_b[100], mA800_g[100], mA800_be[100], mA800_ge[100];

void printArrays(int nbins, bool showErrors=false) {
  if (!showErrors)
    printf("Bin|           other             tt           wlnu        Zleplep        Ztautau          QCDrC             Total\n");
  //.......  1 |          0.3244         0.0157         1.1595         0.0000         7.0269         0.0000 |          8.5265
  else
    printf("Bin|       other                  tt                wlnu             Zleplep             Ztautau               QCDrC                  Total\n");
  //........ 0 |      0.7981+- 0.3566     8.2270+- 1.1635     8.7470+- 3.1302     0.0000+- 0.0000     1.7593+- 0.6514    -2.8657+- 2.8695 |     16.6657+-

  for (int i=0; i<=nbins+1;++i) {
    if (!showErrors) {
      printf("%2i |  %14.4f %14.4f %14.4f %14.4f %14.4f %14.4f |  %14.4f\n", i, 
	     other_sr[i], tt_sr[i], wlnu_sr[i], zleplep_sr[i], ztautau_sr[i], qcd_sr[i],
	     other_sr[i]+ tt_sr[i]+ wlnu_sr[i]+ zleplep_sr[i]+ ztautau_sr[i]+ qcd_sr[i]);
    } else {
      printf("%2i |  %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f |  %10.4f+-%7.4f\n", i, 
	     other_sr[i], other_sr_err[i], tt_sr[i], tt_sr_err[i], wlnu_sr[i], wlnu_sr_err[i],
	     zleplep_sr[i], zleplep_sr_err[i], ztautau_sr[i], ztautau_sr_err[i], qcd_sr[i], qcd_sr_err[i],
	     other_sr[i]+ tt_sr[i]+ wlnu_sr[i]+ zleplep_sr[i]+ ztautau_sr[i]+ qcd_sr[i],
	     sqrt(other_sr_err[i]*other_sr_err[i]+ tt_sr_err[i]*tt_sr_err[i]+ wlnu_sr_err[i]*wlnu_sr_err[i]+ zleplep_sr_err[i]*zleplep_sr_err[i]
		  + ztautau_sr_err[i]*ztautau_sr_err[i]+ qcd_sr_err[i]*qcd_sr_err[i]));
    }
  }
}

void printArraysSignal(int nbins, bool showErrors=false) {
  printf("\n");
  if (!showErrors)
    printf("Bin|BBAprod      200            300            400            500            600            700               800\n");
  //.......  1 |          0.3244         0.0157         1.1595         0.0000         7.0269         0.0000            8.5265
  else
//  printf("Bin|       other                  tt                wlnu             Zleplep             Ztautau               QCDrC                  Total\n");
    printf("Bin|BBAprod  200                 300                 400                 500                 600                 700                    800\n");
  //........ 0 |      0.7981+- 0.3566     8.2270+- 1.1635     8.7470+- 3.1302     0.0000+- 0.0000     1.7593+- 0.6514    -2.8657+- 2.8695 |     16.6657+-
  double sum_200=0, sum_300=0, sum_400=0, sum_500=0, sum_600=0, sum_700=0, sum_800=0;
  for (int i=0; i<=nbins+1;++i) {
    //    if (!showErrors) {
      printf("%2i |  %14.4f %14.4f %14.4f %14.4f %14.4f %14.4f    %14.4f\n", i, 
	     mA200_b[i], mA300_b[i], mA400_b[i], mA500_b[i], mA600_b[i], mA700_b[i], mA800_b[i]);
      sum_200 += mA200_b[i];      sum_300 += mA300_b[i];      sum_400 += mA400_b[i];      sum_500 += mA500_b[i];
      sum_600 += mA600_b[i];      sum_700 += mA700_b[i];      sum_800 += mA800_b[i];
    //    } 
//     else {
//       printf("%2i |  %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f |  %10.4f+-%7.4f\n", i, 
// 	     other_sr[i], other_sr_err[i], tt_sr[i], tt_sr_err[i], wlnu_sr[i], wlnu_sr_err[i],
// 	     zleplep_sr[i], zleplep_sr_err[i], ztautau_sr[i], ztautau_sr_err[i], qcd_sr[i], qcd_sr_err[i],
// 	     other_sr[i]+ tt_sr[i]+ wlnu_sr[i]+ zleplep_sr[i]+ ztautau_sr[i]+ qcd_sr[i],
// 	     sqrt(other_sr_err[i]*other_sr_err[i]+ tt_sr_err[i]*tt_sr_err[i]+ wlnu_sr_err[i]*wlnu_sr_err[i]+ zleplep_sr_err[i]*zleplep_sr_err[i]
// 		  + ztautau_sr_err[i]*ztautau_sr_err[i]+ qcd_sr_err[i]*qcd_sr_err[i]));
//     }
  }
  // plot an extra line with the total sume  
  printf("SUM|  %14.4f %14.4f %14.4f %14.4f %14.4f %14.4f    %14.4f\n", 
	 sum_200, sum_300, sum_400, sum_500, sum_600, sum_700, sum_800);

 sum_200=0, sum_300=0, sum_400=0, sum_500=0, sum_600=0, sum_700=0, sum_800=0;
 printf("\n");
 printf("Bin|GGqAprod      200            300            400            500            600            700               800\n");
  for (int i=0; i<=nbins+1;++i) {
    //    if (!showErrors) {
      printf("%2i |  %14.4f %14.4f %14.4f %14.4f %14.4f %14.4f    %14.4f\n", i, 
	     mA200_g[i], mA300_g[i], mA400_g[i], mA500_g[i], mA600_g[i], mA700_g[i], mA800_g[i]);
      sum_200 += mA200_g[i];      sum_300 += mA300_g[i];      sum_400 += mA400_g[i];      sum_500 += mA500_g[i];
      sum_600 += mA600_g[i];      sum_700 += mA700_g[i];      sum_800 += mA800_g[i];
    //    } 
//     else {
//       printf("%2i |  %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f %10.4f+-%7.4f |  %10.4f+-%7.4f\n", i, 
// 	     other_sr[i], other_sr_err[i], tt_sr[i], tt_sr_err[i], wlnu_sr[i], wlnu_sr_err[i],
// 	     zleplep_sr[i], zleplep_sr_err[i], ztautau_sr[i], ztautau_sr_err[i], qcd_sr[i], qcd_sr_err[i],
// 	     other_sr[i]+ tt_sr[i]+ wlnu_sr[i]+ zleplep_sr[i]+ ztautau_sr[i]+ qcd_sr[i],
// 	     sqrt(other_sr_err[i]*other_sr_err[i]+ tt_sr_err[i]*tt_sr_err[i]+ wlnu_sr_err[i]*wlnu_sr_err[i]+ zleplep_sr_err[i]*zleplep_sr_err[i]
// 		  + ztautau_sr_err[i]*ztautau_sr_err[i]+ qcd_sr_err[i]*qcd_sr_err[i]));
//     }
  }
  // plot an extra line with the total sume  
  printf("SUM|  %14.4f %14.4f %14.4f %14.4f %14.4f %14.4f    %14.4f\n", 
	 sum_200, sum_300, sum_400, sum_500, sum_600, sum_700, sum_800);

}

bool fillArrays(TString name, int i, double cont, double err) {
  if (name =="ATLAS_Bkg_LH_Other") {
    other_sr[i] = cont; other_sr_err[i] = err; return true;
  }
  else if (name == "ATLAS_Bkg_LH_TT") {
    tt_sr[i] = cont; tt_sr_err[i] = err; return true;
  }
  else if (name == "ATLAS_Bkg_LH_Wlnu") {
    wlnu_sr[i] = cont; wlnu_sr_err[i] = err; return true;
  }
  else if (name == "ATLAS_Bkg_LH_Zleplep") {
    zleplep_sr[i] = cont; zleplep_sr_err[i] = err; return true;
  }
  else if (name == "ATLAS_Bkg_LH_Ztautau") {
    ztautau_sr[i] = cont; ztautau_sr_err[i] = err; return true;
  }
  else if (name == "ATLAS_Bkg_LH_qcd_cC") {
    qcd_sr[i] = cont; qcd_sr_err[i] = err; return true;
  }
  //
  else if (name =="ATLAS_Bkg_LH_Other_cB") {
    other_b[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_TT_cB") {
    tt_b[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Wlnu_cB") {
    wlnu_b[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Zleplep_cB") {
    zleplep_b[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Ztautau_cB") {
    ztautau_b[i] = cont;
  }
  //
  else if (name =="ATLAS_Bkg_LH_Other_cD") {
    other_d[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_TT_cD") {
    tt_d[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Wlnu_cD") {
    wlnu_d[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Zleplep_cD") {
    zleplep_d[i] = cont;
  }
  else if (name == "ATLAS_Bkg_LH_Ztautau_cD") {
    ztautau_d[i] = cont;
  }

  // some signal samples
  else if (name == "ATLAS_Sig_LH_bbA100") { mA100_b[i]=cont;  mA100_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA200") { mA200_b[i]=cont;  mA200_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA300") { mA300_b[i]=cont;  mA300_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA400") { mA400_b[i]=cont;  mA400_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA500") { mA500_b[i]=cont;  mA500_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA600") { mA600_b[i]=cont;  mA600_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA700") { mA700_b[i]=cont;  mA700_be[i]=err;}
  else if (name == "ATLAS_Sig_LH_bbA800") { mA800_b[i]=cont;  mA800_be[i]=err;}
  //							     							    
  else if (name == "ATLAS_Sig_LH_ggA100") { mA100_g[i]=cont;  mA100_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA200") { mA200_g[i]=cont;  mA200_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA300") { mA300_g[i]=cont;  mA300_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA400") { mA400_g[i]=cont;  mA400_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA500") { mA500_g[i]=cont;  mA500_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA600") { mA600_g[i]=cont;  mA600_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA700") { mA700_g[i]=cont;  mA700_ge[i]=err;}
  else if (name == "ATLAS_Sig_LH_ggA800") { mA800_g[i]=cont;  mA800_ge[i]=err;}

  return false;
}
/*
ATLAS_Bkg_LH_Other 
ATLAS_Bkg_LH_TT 
ATLAS_Bkg_LH_Wlnu 
ATLAS_Bkg_LH_Zleplep 
ATLAS_Bkg_LH_Ztautau 
ATLAS_Bkg_LH_qcd_cC 

ATLAS_Bkg_LH_Other_cD 
ATLAS_Bkg_LH_TT_cD 
ATLAS_Bkg_LH_Wlnu_cD 
ATLAS_Bkg_LH_Zleplep_cD 
ATLAS_Bkg_LH_Ztautau_cD 

ATLAS_Bkg_LH_Other_cB 
ATLAS_Bkg_LH_TT_cB 
ATLAS_Bkg_LH_Wlnu_cB 
ATLAS_Bkg_LH_Zleplep_cB 
ATLAS_Bkg_LH_Ztautau_cB
*/

// For each shape systematic try to see what is the difference with the original template
void KolmogorovTest(TString fname, double level) {
  
  TFile *f = new TFile(fname, "read");
  TIter next(f->GetListOfKeys());
  TKey *key;
  int nbins_sr = -1;

  vector<TString> v_ks_names_up;
  vector<TString> v_ks_names_do;
  vector<double> v_ks_up;
  vector<double> v_ks_do;




  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1 *h = (TH1*)key->ReadObj();
    TString hname = TString(h->GetName());
    if (hname.Contains("ATLAS_data_")) continue;
    
    if (! hname.Contains("ATLAS_")) continue;
    if (hname.Contains("_cB")) continue;
    if (hname.Contains("_cC") && ! hname.Contains("_qcd_cC")) continue;
    if (hname.Contains("_cD")) continue;

    const int N_syst = 10;
    TString syst[N_syst] = {"_TES", "_JES_EtaModelling", "_JES_FlavComp", "_JES_FlavResp", "_LH_WBKG", 
			    "_MFS", "_TAU_EL",            "_EMBISO",      
			    "_JES_FlavCompG", "_JES_Modelling1"
    };
// JES_EtaModellingup
// JES_FlavCompup
// JES_FlavRespup
// LH_WBKGup
// MFSup
// TAU_ELup
// TESup

    for (int ii=0; ii<N_syst; ++ii) {
      
      if (hname.Contains(syst[ii]+"up")) {
	TString hname_base = hname; hname_base.ReplaceAll(syst[ii]+"up", "");
	TString hname_down = hname; hname_down.ReplaceAll(syst[ii]+"up", syst[ii]+"do");
	if (hname.Contains("MFSup")) {
	  hname_base = "ATLAS_Bkg_LH_Ztautau";
	  hname_down = "ATLAS_Bkg_LH_MFSdo";
	}
	TH1 *h_base = (TH1*) f->Get(hname_base);
	TH1 *h_down = (TH1*) f->Get(hname_down);
	if (!h_base) {cout << "filename: " << hname_base << " not found "  <<  endl;  continue;}
	if (!h_down) {cout << "filename: " << hname_down << " not found "  <<  endl;  continue;}
	if (h->Integral() > 0 && h_base->Integral() > 0 && h_down->Integral() > 0 ) {
	  double ks_up = h_base->KolmogorovTest(h);
	  double ks_do = h_base->KolmogorovTest(h_down);
	  printf("# Histo %25s TES  KSup= %6.4f  KSdo= %6.4f\n", hname_base.Data(), ks_up, ks_do);
	  if (ks_up>level && ks_do>level) {
	    v_ks_names_up.push_back(hname); v_ks_names_do.push_back(hname_down);
	    v_ks_up.push_back(ks_up); v_ks_do.push_back(ks_do);
	  }
	}
      }

    }
    
  }
  cout << endl;
  cout << "# List of Histos with both up and down level > " << level << " for file: " << fname  << endl;

  for (int i=0 ;i < v_ks_names_up.size(); ++i) {
    if (  (v_ks_names_up.at(i)).Contains("_Sig_")  ) continue;
    printf(">>> %15s # %35s  up= %6.4f do= %6.4f\n", fname.Data(), (v_ks_names_up.at(i)).Data(), v_ks_up.at(i), v_ks_do.at(i));
  }
  cout << endl;
  for (int i=0 ;i < v_ks_names_up.size(); ++i) {
    if (  (v_ks_names_up.at(i)).Contains("_Sig_") || (v_ks_names_up.at(i)).Contains("_XSec_")  ) 
      printf(">>> %15s # %35s  up= %6.4f do= %6.4f\n", fname.Data(), (v_ks_names_up.at(i)).Data(), v_ks_up.at(i), v_ks_do.at(i));
  }
  cout << endl;
  //  cout << "sed  ";  
  for (int i=0 ;i < v_ks_names_up.size(); ++i) {
    cout <<  "sed  -i \"/"<< v_ks_names_up.at(i)  << "/d\" $FILES" << endl;
    cout <<  "sed -i \"/"<< v_ks_names_do.at(i)  << "/d\" $FILES " << endl;
  }
  cout << endl; 
}


// compare integrals for 2 files with the same input
void FileComparison(TString fname1, TString fname2) {
  TFile *f1 = new TFile(fname1, "read");
  TFile *f2 = new TFile(fname2, "read");

  std::map<TString,double> ints_f1;
  std::map<TString,double> ints_f2;

  TIter next(f1->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1 *h = (TH1*)key->ReadObj();
    TString hname = TString(h->GetName());
    if (hname.Contains("ATLAS_data_")) continue; 
    ints_f1[hname] = h->Integral();
  }
  
  next= f2->GetListOfKeys();
  TKey *key2;
  while ((key2 = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key2->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1 *h = (TH1*)key2->ReadObj();
    TString hname = TString(h->GetName());
    if (hname.Contains("ATLAS_data_")) continue; 
    ints_f2[hname] = h->Integral();
  }
  
  // now print the signal
  int const n_sig_mass = 11;
  TString sig_mass[n_sig_mass]={"090","100","110","120","130","140","150","170","200","300","500"};

  for (int i=0; i < n_sig_mass; ++i) {
    TString name = "ATLAS_Sig_LH_bbA" + sig_mass[i];
    double c_f1 = ints_f1.find(name)->second;
    double c_f2 = ints_f2.find(name)->second;
    TString name1 = "ATLAS_Sig_LH_ggA" + sig_mass[i];
    double c_f3 = ints_f1.find(name1)->second;
    double c_f4 = ints_f2.find(name1)->second;

    printf("bbA%3s  %10.4f  %10.4f  %10.2f     ggA%3s %10.4f  %10.4f  %10.2f\n", sig_mass[i].Data(), 
	   c_f1, c_f2, c_f1/c_f2, sig_mass[i].Data(), c_f3, c_f4, c_f3/c_f4);

  }
  
  // same for backgrounds
  int const n_bkgs = 6;
  TString bkgs[n_bkgs]={"qcd_cC","Ztautau","Wlnu","Zleplep","Other","TT"};
  
  for (int i=0; i < n_bkgs; ++i) {
    TString name = "ATLAS_Bkg_LH_" + bkgs[i];
    double c_f1 = ints_f1.find(name)->second;
    double c_f2 = ints_f2.find(name)->second;
    if (i == 0) {
      printf("%7s  %10.4f  %10.4f  %10.2f\n", bkgs[i].Data(), c_f1, c_f2, c_f1/c_f2);
    } else {
      double c_f1_b = ints_f1.find(name+"_cB")->second;
      double c_f2_b = ints_f2.find(name+"_cB")->second;
      double c_f1_c = ints_f1.find(name+"_cC")->second;
      double c_f2_c = ints_f2.find(name+"_cC")->second;
      double c_f1_d = ints_f1.find(name+"_cD")->second;
      double c_f2_d = ints_f2.find(name+"_cD")->second;
      
      printf("%7s  %10.4f  %10.4f  %10.2f | %10.4f  %10.4f  %10.2f | %10.4f  %10.4f  %10.2f \
| %10.4f  %10.4f  %10.2f \n", 
	     bkgs[i].Data(), c_f1, c_f2, c_f1/c_f2,  
	     c_f1_b, c_f2_b, c_f1_b/c_f2_b,  c_f1_c, c_f2_c, c_f1_c/c_f2_c,  c_f1_d, c_f2_d, c_f1_d/c_f2_d);
      
    }

  }
}

void HistoReader(TString fname, bool showerrors=false, bool plotSignal=false) {

  TFile *f = new TFile(fname, "read");

//   std::map<int, double> BinSum_CR_B;
//   std::map<int, double> BinSum_CR_D;
//   bool mapCreated = false;
  TIter next(f->GetListOfKeys());
  TKey *key;
  int nbins_sr = -1;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1 *h = (TH1*)key->ReadObj();
    //    cout << h->GetName() << endl;
    TString hname = TString(h->GetName());
    //bool isSignal = hname.Contains("ATLAS_Sig") || hname.Contains("ggA") || hname.Contains("bbA");
    if (hname.Contains("ATLAS_data_")) continue;
    //if (hname.EndsWith("up")) continue;
    //if (hname.EndsWith("do")) continue;
    //cout << "test" << endl;
    int nbins = h->GetNbinsX();
//     if (!mapCreated) {
//       for (int i=0; i<nbins+1; ++i) {
// 	BinSum[i] = 0.;
//       }
//       mapCreated = true;
//     }
    if (h->Integral() <=0) {
      cout << "ERROR: histo " << hname << " Integral = " << h->Integral() << endl;
    }
    for (int i=0; i <= nbins+1; ++i) {
      double h_bin = h->GetBinContent(i);
      double h_err = h->GetBinError(i);
      if (fillArrays(hname, i, h_bin, h_err)) nbins_sr = nbins;
      if (i==0 && h_bin != 0) {
	//cout << "histo: " << h->GetName() << " has non zero underflow: " << h_bin  << endl;
      }
      else if (i==nbins+1 && h_bin != 0) {
	//cout << "histo: " << h->GetName() << " has non zero overflow: " << h_bin << endl;
      }
      else if (h_bin < 0) {
	cout << "Bin " << i  << " histo: " << h->GetName() <<  " entries: " << h_bin << endl;
      }
    }

    
  }
  printArrays(nbins_sr,showerrors);
  if (plotSignal) {
    printArraysSignal(nbins_sr,showerrors);
  }
}

TH1F* getModifiedHisto(TH1F* h, bool USEKYLE=false) {
  int nbins = h->GetNbinsX();
  TString hname= h->GetName();
  TString htitle= h->GetTitle();
  double xmin=h->GetBinLowEdge(1);
  double xmax=h->GetBinLowEdge(nbins)+h->GetBinWidth(nbins);

  // try to get ingredients needed for Kyle's prescription
  // average weight
  double averageErr2 = 0;
  for (int i=1; i<= h->GetEntries(); ++i) {
    averageErr2 += h->GetBinError(i)*h->GetBinError(i);
  }
  double averageWeightError = sqrt(averageErr2/h->GetEntries());
  double averageWeight = h->GetSumOfWeights()/h->GetEntries();

  TH1F * hh = new TH1F(hname, htitle, nbins, xmin, xmax);

  hh->SetBinContent( 0, h->GetBinContent(0) );
  hh->SetBinError( 0, h->GetBinError(0) );

  hh->SetBinContent( nbins+1, h->GetBinContent(nbins+1) );
  hh->SetBinError( nbins+1, h->GetBinError(nbins+1) );
  cout << "Working on histo " << hname << endl;
  for (int i=1; i <= nbins; ++i) {
    double h_bin = h->GetBinContent(i);
    double h_err = h->GetBinError(i);
    if (h_bin > 0) {
      hh->SetBinContent(i, h_bin);
      hh->SetBinError(i, h_err);      
    }
    else if (h_bin==0 && USEKYLE) { // if this is 0 but not negative, implement Kyle's prescription
      hh->SetBinContent(i, averageWeight);
      hh->SetBinError(i, averageWeightError);      
    }
    else {
      // this will probably don't work - try to use Kyle's prescription
//       hh->SetBinContent(i, 0);
//       hh->SetBinError(i, sqrt(h_bin*h_bin + h_err*h_err));
      hh->SetBinContent(i, averageWeight);
      hh->SetBinError(i, averageWeightError);      

    }
    printf("Average Weights: %12.8f +- %12.8f\n",averageWeight, averageWeightError);
    printf("%2i  %14.4f +- %8.4f   %14.4f +- %8.4f \n", i, h_bin, h_err, hh->GetBinContent(i), hh->GetBinError(i));
  }
  return hh;
}

// ---- code for systematics smoothing ---------- //
TH1F* EqualAreaGabriel(TH1F* n, TH1F* s) {
  //cout << "EQUAL AREA" << endl;
  float frac = 0.05; // minimal statistical error

  // find boundaries where integral = frac of max
  vector< pair<int, int> > bins;

  // find max of first set
  int maxBin(23);
  float tot(0);
  float err(0);

  for(int b=1; b<n->GetNbinsX()+1; b++) {
    //cout << b <<  "\t" << maxBin << "\t" << n->GetBinError(b) / n->GetBinContent(b) << endl;
    if(n->GetBinError(b) / n->GetBinContent(b) > frac){
      tot = 0;
      err = 0;
      for(int c=b; c<maxBin+1; c++) {
        tot += n->GetBinContent(c);
        err = sqrt(err*err + n->GetBinError(c)*n->GetBinError(c));
        //cout << "\t" << err / tot << endl;
        // expanding integral until cover fraction & picks up tails
        if( err / tot < frac ) {
          // need to catch tails of distribution and do not want to cut off too soon
          // if integral of remaining bins is < frac of max && that + region in question < max - take it all
          float tmpTot(0);
          float tmpErr(0);
          for(int d=c+1; d<maxBin+1; d++) {
            tmpTot += n->GetBinContent(d);
            tmpErr = sqrt(tmpErr*tmpErr + n->GetBinError(d)*n->GetBinError(d));
          }
          if( tmpErr / tmpTot > frac) {
            c = maxBin;
          }

          bins.push_back( make_pair(b,c) );
          b = c; // advance outer loop
          break;
        }

        // save last pair before loop completes
        if(c==maxBin){
          bins.push_back( make_pair(b,c) );
          b=c;
          break;
        } // if loop ending
      } // second loop
    } // if content < frac of max
    else{
      bins.push_back( make_pair(b,b) );
    }
    if(b==maxBin) { // in the last bin
      maxBin += 23;
    }
  } // outer loop over bins

  //  TString newName("new_");
  //  newName.Append(s->GetName());
  TString newName(s->GetName());
  TH1F* newSyst = new TH1F(newName, newName, s->GetNbinsX(), s->GetBinLowEdge(1), s->GetBinLowEdge( s->GetNbinsX()+1 ));
  newSyst->Sumw2();

  for(unsigned int i=0; i<bins.size(); i++) {
    float nom=n->Integral( (bins.at(i)).first, (bins.at(i)).second );
    float sys=s->Integral( (bins.at(i)).first, (bins.at(i)).second );
    // fill new systematic hist
    for(int b=(bins.at(i)).first; b<(bins.at(i)).second+1; b++) {
      if(nom>0) {
        newSyst->SetBinContent(b, (sys/nom) * n->GetBinContent(b));
      } else {
        newSyst->SetBinContent(b, nom);
      }
    }
  } // loop over equal area bins
  return newSyst;
} // EqualArea




// ---------------------------------------------- //
// end: systematics smoothing code

// read this file and update it by addind a histogram
// 
void AddEMBISODown(TString infile) {
  
  TFile *f = new TFile(infile, "update");

  TH1F *h    = (TH1F*) f->Get("ATLAS_Bkg_LH_Ztautau");
  TH1F *h_up = (TH1F*) f->Get("ATLAS_Bkg_LH_Ztautau_EMBISOup");

  TH1F *h_do = new TH1F("ATLAS_Bkg_LH_Ztautau_EMBISOdo_SYMM","ATLAS_Bkg_LH_Ztautau_EMBISOdo_SYMM", 
			h->GetNbinsX(), h->GetBinLowEdge(1), 
			h->GetBinLowEdge( h->GetNbinsX() ) + h->GetBinWidth(h->GetNbinsX()));
  cout << "Embedding Isolation Down syst:  file " << infile << endl;

  for (int i=1; i <= h->GetNbinsX(); ++i) {
    double x_no = h->GetBinContent(i);
    double x_up = h_up->GetBinContent(i);
    double rel_var = (x_no!=0)?(x_up - x_no) / x_no : 0.;
    double x_do = x_no * ( 1 - rel_var );
    h_do->SetBinContent(i,x_do);
    printf("Nominal: %10.3f   Up: %10.3f / %7.4f    Down: %10.3f \n", x_no, x_up, rel_var, x_do);
  }
  
  h_do->Write("ATLAS_Bkg_LH_Ztautau_EMBISOdo_SYMM");
  f->Close();
  
}

void UpdateHistos(TString infile, TString outfile, bool USEKYLE=false) {
  TH1F::SetDefaultSumw2();
  TFile *f = new TFile(infile, "read");
  if (outfile == "") outfile = "upd_"+infile;
  TFile *fout = new TFile(outfile, "recreate");
  
  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    int nbins = h->GetNbinsX();
    // check if there is a bin with negative entries
    bool hasNegativeEntries=false;
    for (int i=1; i <= nbins; ++i) {
      double h_bin = h->GetBinContent(i);
      bool hasneg = USEKYLE ? h_bin <= 0 : h_bin < 0;
      if (hasneg) { hasNegativeEntries = true; break;  } // include also zero entries for Kyle's prescription
    }
    if (hasNegativeEntries) {
      TH1F *hh = getModifiedHisto(h, USEKYLE);      
      hh->Write();
    } else {
      h->Write();
    }

  }
  fout->Close();
}

void SmoothSystematics(TString infile, TString outfile) {

  // for this to work we need to know the names of the nominal histograms
  const int n_bkgList = 6+42;
  TString bkgList[n_bkgList] ={"ATLAS_Bkg_LH_Other", "ATLAS_Bkg_LH_TT", "ATLAS_Bkg_LH_Wlnu",
			       "ATLAS_Bkg_LH_Zleplep", "ATLAS_Bkg_LH_Ztautau", "ATLAS_Bkg_LH_qcd_cC",
			       "ATLAS_Sig_LH_bbA090",  "ATLAS_Sig_LH_bbA090",
			       "ATLAS_Sig_LH_bbA100",  "ATLAS_Sig_LH_bbA100",
			       "ATLAS_Sig_LH_bbA110",  "ATLAS_Sig_LH_bbA110",
			       "ATLAS_Sig_LH_bbA120",  "ATLAS_Sig_LH_bbA120",
			       "ATLAS_Sig_LH_bbA125",  "ATLAS_Sig_LH_bbA125",
			       "ATLAS_Sig_LH_bbA130",  "ATLAS_Sig_LH_bbA130",
			       "ATLAS_Sig_LH_bbA140",  "ATLAS_Sig_LH_bbA140",
			       "ATLAS_Sig_LH_bbA150",  "ATLAS_Sig_LH_bbA150",
			       "ATLAS_Sig_LH_bbA170",  "ATLAS_Sig_LH_bbA170",
			       "ATLAS_Sig_LH_bbA200",  "ATLAS_Sig_LH_bbA200",
			       "ATLAS_Sig_LH_bbA250",  "ATLAS_Sig_LH_bbA250",
			       "ATLAS_Sig_LH_bbA300",  "ATLAS_Sig_LH_bbA300",
			       "ATLAS_Sig_LH_bbA350",  "ATLAS_Sig_LH_bbA350",
			       "ATLAS_Sig_LH_bbA400",  "ATLAS_Sig_LH_bbA400",
			       "ATLAS_Sig_LH_bbA450",  "ATLAS_Sig_LH_bbA450",
			       "ATLAS_Sig_LH_bbA500",  "ATLAS_Sig_LH_bbA500",
			       "ATLAS_Sig_LH_bbA550",  "ATLAS_Sig_LH_bbA550",
			       "ATLAS_Sig_LH_bbA600",  "ATLAS_Sig_LH_bbA600",
			       "ATLAS_Sig_LH_bbA700",  "ATLAS_Sig_LH_bbA700",
			       "ATLAS_Sig_LH_bbA800",  "ATLAS_Sig_LH_bbA800",
			       "ATLAS_Sig_LH_bbA900",  "ATLAS_Sig_LH_bbA900"

};
  //  "ATLAS_Sig_LH_bbA150"

  TH1F::SetDefaultSumw2();
  TFile *f = new TFile(infile, "read");
  if (outfile == "") outfile = "upd_"+infile;
  TFile *fout = new TFile(outfile, "recreate");
  
  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    // check whether this is an upwards variation
    bool thisIsNotVariation=true;
    for (int i=0; i < n_bkgList; ++i) {
      if (hname.Contains(bkgList[i]) && (hname.EndsWith("up") || hname.EndsWith("do"))) {
	// get the nominal histo
	TH1F *h_nominal = (TH1F*) f->Get(bkgList[i]);
	TH1F *hh = EqualAreaGabriel(h_nominal, h);
	hh->Write(hname);

	thisIsNotVariation=false;
	cout << "Found variation histogram " << hname << " and comparing the changes: " << endl;
	for (int jj=1; jj <= h->GetNbinsX(); ++jj) {
	  if (h->GetBinContent(jj) > 0) {
	    printf("bin %2i  %8.3f   %8.3f   %8.3f   %8.3f  %8.3f %% \n", jj,
		   h_nominal->GetBinContent(jj), h->GetBinContent(jj), hh->GetBinContent(jj),
		   h->GetBinContent(jj) - hh->GetBinContent(jj),  
		   (h->GetBinContent(jj) - hh->GetBinContent(jj))*100./h->GetBinContent(jj));

	  } else {
	    printf("bin %2i  %8.3f   %8.3f   %8.3f   %8.3f \n", jj,
		   h_nominal->GetBinContent(jj), h->GetBinContent(jj), hh->GetBinContent(jj),
		   h->GetBinContent(jj) - hh->GetBinContent(jj) );
	  }
	}

	break;
      }
    }
    // if this is not a variation histogram, simply write what you have gotten!
    if (thisIsNotVariation) {
      h->Write();
    }


  }
  fout->Close();

}

// function: insert the name of a NP and check for all shape systematics
//           what is the KS test value for the nominal
//           also make the plots for the + 1 sigma and - 1 sigma variations
//           of this NP in the total sum of bkgs (no fitting involved)
void CheckKS(TString infile, TString inNP) {
  
  TFile *f = new TFile(infile, "read");
  TIter next(f->GetListOfKeys());
  TKey *key;

  cout << "Working for file: " << infile << endl;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    if (! hname.Contains("ATLAS_Bkg")) continue;
    if ( (hname.Contains("_cC") && !hname.Contains("qcd")) || hname.Contains("_cD") 
	 || hname.Contains("_cB") || hname.Contains("_data_") || hname.Contains("AlpZtt") ) continue;

    if (hname.Contains( (inNP+"up") )) {
      // cout << "Found: " << hname << endl;
      TString hname_down = hname;
      hname_down.ReplaceAll( (inNP+"up"), (inNP+"do"));
      TString hname_nominal = hname;
      hname_nominal.ReplaceAll( ("_"+inNP+"up"), "");
      // cout << "Nominal: " << hname_nominal << endl;
      // cout << "Down:    " << hname_down << endl;

      TH1F *h_do = (TH1F*) f->Get(hname_down);
      TH1F *h_no = (TH1F*) f->Get(hname_nominal);
      //cout << "Check nominal: " << h_no << endl;
      double KS_up = h_no->KolmogorovTest(h, "M");
      double KS_do = h_no->KolmogorovTest(h_do, "M");
      double KS_up_prob = h_no->KolmogorovTest(h);
      double KS_do_prob = h_no->KolmogorovTest(h_do);
      // double chi2_up = h_no->Chi2Test(h, "WW");
      // double chi2_do = h_no->Chi2Test(h_do, "WW");

      // printf("%20s %10s  up: %10.8f / %10.8f  do: %10.8f / %10.8f \n", 
      // 	     hname_nominal.Data(), inNP.Data(), KS_up, chi2_up, KS_do, chi2_do);
      printf("%20s %10s  up: %10.8f (%10.8f)  do: %10.8f (%10.8f)  %s %s\n", 
	     hname_nominal.Data(), inNP.Data(), KS_up, KS_up_prob,
	     KS_do, KS_do_prob, hname.Data(), hname_down.Data());

      // try to test the differences bin by bin
      int nbins = h_no->GetNbinsX();
      TGraph * g_nominal_cumulative = new TGraph();
      g_nominal_cumulative->SetPoint(0, 0, 0);
      double nominal_integral =  h_no->Integral(0,999);
      for (int i=1; i<=nbins; ++i) {
	g_nominal_cumulative->SetPoint(i, i, h_no->Integral(0,i)/nominal_integral);
      }
      double D_max = -999;
      for (int i=1; i<=nbins; ++i) {
	double cont_nom = h_no->GetBinContent(i);
	double diff_up = (cont_nom>0)?(h->GetBinContent(i) - cont_nom)/cont_nom:0.;
	double diff_do = (cont_nom>0)?(h_do->GetBinContent(i) - cont_nom)/cont_nom:0.;
	//	printf("%2i up: %10.7f  do: %10.7f \n",  i, diff_up, diff_do);
	// try to show the cumulative distribution
	double F_i_do = g_nominal_cumulative->Eval( h_do->Integral(0,i)/nominal_integral  );
	double F_i_up = g_nominal_cumulative->Eval( h->Integral(0,i)/nominal_integral  );
	double D_i_up = fabs(h_no->Integral(0,i)/nominal_integral -h->Integral(0,i)/nominal_integral ) ;
	double D_i_up_x = std::max(F_i_up-(i-1)/nbins, i/nbins - F_i_up);
	printf("%2i up: %10.5f  max: %10.5f %10.5f \n",  i, cont_nom, D_i_up,  D_i_up_x);
	if (D_i_up > D_max) D_max = D_i_up;
      }
      cout << "D_max = " << D_max << endl;
    }

  }
  

  
  
  
}



// Test function: rebin all the histograms in a particular file to a user defined scheme
void fileRebin(TString infile, int nrebin) {
  
  TFile *f = new TFile(infile, "read");
  TString outfile = "rebinned_";
  outfile+=nrebin;
  outfile+=("_"+infile);
  TFile *fout = new TFile(outfile, "recreate");

  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    h->Rebin(nrebin);
    h->Write(hname);
  }
  
  fout->Close();
}
// 
// special rebinning for 31 bins in the histo
// 
void fileRebinSpecial(TString infile, int nrebin) {
  
  TFile *f = new TFile(infile, "read");
  TString outfile = "rebinned_";
  outfile+=nrebin;
  outfile+=("_"+infile);
  TFile *fout = new TFile(outfile, "recreate");

  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    TH1F *hnew = new TH1F(); hnew->SetName(hname);
    if (nrebin == 4) {
      double x_rebin[5] = {0,10, 15,20,30};
      h->Rebin(nrebin, "hnew", x_rebin);
    }
    else    if (nrebin == 3000) {
      double x_rebin[4] = {0,10,20,31};
      h->Rebin(nrebin, "hnew", x_rebin);
    }
    else  if (nrebin == 3) {
      double x_rebin[4] = {0,15,31,50};
      h->Rebin(nrebin, "hnew", x_rebin);
    }
    // else {
    //   h->Rebin(nrebin);
    // }
    // add the overflow into the last bin and set the overflow to zero
    hnew->SetName(hname);
    hnew->Write(hname);
  }
  
  fout->Close();
}

// ===========================================================================================
//
// extra function: replace each signal histogram with a gaussian with the same
// mean and RMS to test the shape effects
//
TH1F* getRealBinHisto(TH1F* h, TH1F *lowEdge) {
  double *binning = new double[lowEdge->GetNbinsX()+1];
  
  binning[0] = 0;
  for (int i=1; i <lowEdge->GetNbinsX(); ++i){
    binning[i] = lowEdge->GetBinContent(i);
  }
  binning[lowEdge->GetNbinsX()] = lowEdge->GetBinContent( lowEdge->GetNbinsX()-1 )+100;

  TH1F *h1_new = new TH1F(TString(h->GetName())+TString("_real"), TString(h->GetName()), lowEdge->GetNbinsX(), binning);

  for (int i=1; i<=lowEdge->GetNbinsX(); ++i) {
    h1_new->SetBinContent(i, h->GetBinContent(i));
    h1_new->SetBinError(i, h->GetBinError(i));
  }

  delete [] binning;
  return h1_new;
}

TH1F* getFakeBinHisto(TH1F* h) {
  
  TH1F *h1_new = new TH1F(TString(h->GetName())+TString("_fake"), TString(h->GetName()), h->GetNbinsX(), 0, double(h->GetNbinsX()));

  for (int i=1; i<=h->GetNbinsX(); ++i) {
    h1_new->SetBinContent(i, h->GetBinContent(i));
    h1_new->SetBinError(i, h->GetBinError(i));
  }
  return h1_new;
}

TH1F *getGaussianHisto(TH1F *h) {

  double mean = h->GetMean();
  double sigma = h->GetRMS()/2.;
  double integral = h->Integral();

  TH1F *h_new = new TH1F(*h);
  h_new->SetName(TString(h->GetName())+TString("_gauss"));
  for (int i=1; i<= h->GetNbinsX(); ++i) {
    double x = h_new->GetBinCenter(i);
    h_new->SetBinContent(i, TMath::Gaus(x, mean, sigma));
    h_new->SetBinError(i, 0);
  }
  h_new->Scale( integral /h_new->Integral() );
  
  return h_new;

}

TH1F* getModifiedSignalHisto(TH1F* h, TH1F *lowEdge) {
  TString hname= h->GetName();
  TString htitle= h->GetTitle();
  // get the real histo:
  TH1F *h_real = getRealBinHisto(h, lowEdge);
  // get a Gaussian histo based on the real histo
  TH1F *h_real_gaussian = getGaussianHisto(h_real);
  // now convert the real gaussian to a fake one
  TH1F *h_fake_gaussian = getFakeBinHisto(h_real_gaussian);
  cout << "TEST: default " << h->Integral() << endl;
  cout << "TEST: real    " << h_real->Integral() << endl;
  cout << "TEST: realg   " << h_real_gaussian->Integral() << endl;
  cout << "TEST: fakeg   " << h_fake_gaussian->Integral() << endl;

  return h_fake_gaussian;
}

void UpdateHistosAndSignal(TString infile, TString outfile) {

  TH1F::SetDefaultSumw2();
  TFile *f = new TFile(infile, "read");
  if (outfile == "") outfile = "upd_"+infile;
  TFile *fout = new TFile(outfile, "recreate");
  
  TH1F *lowEdge = (TH1F*) f->Get("Bin_Lower_Edges");
  
  TIter next(f->GetListOfKeys());
  TKey *key;
  while ((key = (TKey*)next())) {
    TClass *cl = gROOT->GetClass(key->GetClassName());
    if (!cl->InheritsFrom("TH1")) continue;
    TH1F *h = (TH1F*)key->ReadObj();
    TString hname = TString(h->GetName());
    int nbins = h->GetNbinsX();
    // check if there is a bin with negative entries
    bool hasNegativeEntries=false;
    for (int i=1; i <= nbins; ++i) {
      double h_bin = h->GetBinContent(i);
      if (h_bin < 0) { hasNegativeEntries = true; break;  }
    }
    bool isSignal=false;
    if (hname.Contains("Sig_LH")) { 
      isSignal=true;
    }
    if (isSignal) {
      TH1F *hh = getModifiedSignalHisto(h, lowEdge);      
      hh->Write(h->GetName());
    }
    else if (hasNegativeEntries) {
      TH1F *hh = getModifiedHisto(h);      
      hh->Write();
    } 
    else {
      h->Write();
    }

  }
  fout->Close();
}