#include <TFile.h>
#include <TDirectory.h>
#include <TH2.h>
#include <TAxis.h>
#include <TMath.h>
#include <TStyle.h>
#include <TCanvas.h>

void Add(TH2 *hout, const TH2 *hin) {
  if ((!hout) || (!hin)) return; // just a precaution
  if (!(hout->GetSumw2N() > 0)) hout->Sumw2(kTRUE);
  int flow = 0; // 0 = "include", 1 = "exclude" under/over-flow bins
  for (int i = flow; i <= (hin->GetNbinsX() - flow + 1); i++) {
    double x = hin->GetXaxis()->GetBinCenter(i);
    for (int j = flow; j <= (hin->GetNbinsY() - flow + 1); j++) {
      double c = hin->GetBinContent(i, j);
      double w = hin->GetBinError(i, j);
      if ((c != 0.0) || (w != 0.0)) {
        double y = hin->GetYaxis()->GetBinCenter(j);
        int n = hout->FindFixBin(x, y);
        
        // hout->AddBinContent(n, c);
        hout->SetBinContent(n, (hout->GetBinContent(n) + c));
        hout->SetBinError(n, TMath::Sqrt(c * c + w * w));
      }
    }
  }
  hout->SetEntries(hin->GetEntries());
  hout->ResetStats();
}

void bins_zjets_norm() {
  TFile *f1 = TFile::Open("MC13TeV_DY50toInf_nlo.root");
  TH2 *h1_orig; gDirectory->GetObject("third_tight_mu_pt_eta", h1_orig);
  TH2 *h2_orig; gDirectory->GetObject("third_loose_mu_pt_eta", h2_orig);
  TH2 *nt2_orig; gDirectory->GetObject("nt2_pt_eta", nt2_orig);
  
double scale = 1921.8;

h1_orig->Scale(scale);
h2_orig->Scale(scale);
nt2_orig->Scale(scale);

#if 1 /* 0 or 1 */
  // "coinciding" bin edges
  double bins_x[] = {
    h1_orig->GetXaxis()->GetXmin(),
//    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(10.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(20.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(40.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(70.) + 1),
//    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(100.) + 1),
    h1_orig->GetXaxis()->GetXmax()
  };
  double bins_y[] = {
    h1_orig->GetYaxis()->GetXmin(),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(0.)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(1.)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(1.5)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(2.5) + 1),
    h1_orig->GetYaxis()->GetXmax()
  };
  int nxf = 2, nyf = 2; // the first bin to use / draw
  int nxl = (sizeof(bins_x) / sizeof(double) - 2); // the last bin to use / draw
  int nyl = (sizeof(bins_y) / sizeof(double) - 2); // the last bin to use / draw
#elif 1 /* 0 or 1 */
  // "coinciding" bin edges
  double bins_x[] = {
//    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(10.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(20.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(40.)),
    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(70.) + 1),
//    h1_orig->GetXaxis()->GetBinLowEdge(h1_orig->GetXaxis()->FindFixBin(100.) + 1)
  };
  double bins_y[] = {
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(0.)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(1.)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(1.5)),
    h1_orig->GetYaxis()->GetBinLowEdge(h1_orig->GetYaxis()->FindFixBin(2.5) + 1)
  };
  int nxf = 0, nyf = 0; // use / draw all bins
  int nxl = 0, nyl = 0; // use / draw all bins
#elif 1 /* 0 or 1 */
  // "arbitrary" bin edges
  double bins_x[] = {0., 10., 20., 40., 70., 100., 300.};
  double bins_y[] = {-3., 0., 1., 1.5, 2.5, 3.};
  int nxf = 2, nyf = 2; // the first bin to use / draw
  int nxl = (sizeof(bins_x) / sizeof(double) - 2); // the last bin to use / draw
  int nyl = (sizeof(bins_y) / sizeof(double) - 2); // the last bin to use / draw
#else /* 0 or 1 */
  // "arbitrary" bin edges
  double bins_x[] = {10., 20., 40., 70., 100.};
  double bins_y[] = {0., 1., 1.5, 2.5};
  int nxf = 0, nyf = 0; // use / draw all bins
  int nxl = 0, nyl = 0; // use / draw all bins
#endif /* 0 or 1 */
  
  TH2F *h1 = new TH2F("h1", "tight muons", (sizeof(bins_x) / sizeof(double) - 1), bins_x, (sizeof(bins_y) / sizeof(double) - 1), bins_y);
  h1->Sumw2(kTRUE);
  Add(h1, h1_orig);
  
  TH2F *h2 = new TH2F("h2", "loose muons", (sizeof(bins_x) / sizeof(double) - 1), bins_x, (sizeof(bins_y) / sizeof(double) - 1), bins_y);
  h2->Sumw2(kTRUE);
  Add(h2, h2_orig);
  
  TH2F *nt2 = new TH2F("nt2", "nt2", (sizeof(bins_x) / sizeof(double) - 1), bins_x, (sizeof(bins_y) / sizeof(double) - 1), bins_y);
  nt2->Sumw2(kTRUE);
  Add(nt2, nt2_orig);
  
#if 1 /* 0 or 1 */
  h1->GetXaxis()->SetRange(nxf, nxl);
  h1->GetYaxis()->SetRange(nyf, nyl);
  h2->GetXaxis()->SetRange(nxf, nxl);
  h2->GetYaxis()->SetRange(nyf, nyl);
  nt2->GetXaxis()->SetRange(nxf, nxl);
  nt2->GetYaxis()->SetRange(nyf, nyl);
#endif /* 0 or 1 */
  
  TH2F *f = new TH2F(*h1);
  f->SetNameTitle("f", "Fake Ratio");
  f->Divide(h2);
  
  double mul = 1.0;
#if 1 /* 0 or 1 */
  // use bins in the set bin ranges
  // if (h2->GetEffectiveEntries() != 0.0) mul = h1->GetEffectiveEntries() / h2->GetEffectiveEntries();
  if (h2->Integral() != 0.0) mul = h1->Integral() / h2->Integral();
#else /* 0 or 1 */
  // use all entries
  if (h2->GetEntries() != 0.0) mul = h1->GetEntries() / h2->GetEntries();
#endif /* 0 or 1 */
  if (mul != 1.0) mul /= (1.0 - mul);
  
  TH2F *Nt3 = new TH2F(*nt2);
  Nt3->SetNameTitle("Nt3", "Nt3");
  Nt3->Scale(mul);
  
  gStyle->SetOptStat("neMRuoI");
  TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
//  c1->Divide(2, 2);
//  c1->cd(1);
//  f->Draw("colz text e");
//  c1->cd(2);
//  h1->Draw("colz text e");//tight
//  h2->Draw("colz text e");//loose

//  c1->cd(3);
  nt2->Draw("colz text e");
//  c1->cd(4);
//  Nt3->Draw("colz text e");
//c1->SaveAs("zjets_tight_norm.root");
//c1->SaveAs("zjets_loose_norm.root");
//c1->SaveAs("zjets_fr_norm.root");
c1->SaveAs("zjets_nt2_norm.root");
//c1->SaveAs("zjets_nt3_norm.root");
}