#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 new_bins() {
  TFile *f1 = TFile::Open("MC13TeV_TTJets.root");
  TH2 *h1_orig; f1->GetObject("third_tight_mu_pt_eta", h1_orig);
  TH2 *h2_orig; f1->GetObject("third_loose_mu_pt_eta", h2_orig);
//  TH2 *nt2t_or  TFile *f1 = TFile::Open("MC13TeV_TTJets.root");
  TH2 *nt2_orig; f1->GetObject("nt2_pt_eta", nt2_orig);


  TFile *f2 = TFile::Open("MC13TeV_DY50toInf_nlo.root");
  TH2 *h1d_orig; f2->GetObject("third_tight_mu_pt_eta", h1d_orig);
  TH2 *h2d_orig; f2->GetObject("third_loose_mu_pt_eta", h2d_orig);
  TH2 *nt2d_orig; f2->GetObject("nt2_pt_eta", nt2d_orig);

   TFile *f3 = TFile::Open("MC13TeV_WWTo2L2Nu.root");
  TH2 *h1w_orig; f2->GetObject("third_tight_mu_pt_eta", h1w_orig);
  TH2 *h2w_orig; f2->GetObject("third_loose_mu_pt_eta", h2w_orig);
  TH2 *nt2w_orig; f2->GetObject("nt2_pt_eta", nt2w_orig);
 
double scale1 = 0.731871, scale2 = 0.000107827, scale3 = 0.0243562;

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

h1d_orig->Scale(scale2);
h2d_orig->Scale(scale2);
nt2d_orig->Scale(scale2);

h1w_orig->Scale(scale3);
h2w_orig->Scale(scale3);
nt2w_orig->Scale(scale3);

h1_orig->Add(h1d_orig);
h1_orig->Add(h1w_orig);
h2_orig->Add(h2d_orig);
h2_orig->Add(h2w_orig);
nt2_orig->Add(nt2d_orig);
nt2_orig->Add(nt2w_orig);

#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");

//  c1->cd(3);
//  nt2->Draw("colz text e");
//  c1->cd(4);
//  Nt3->Draw("colz text e");
//c1->SaveAs("added_tight_norm.root");
//c1->SaveAs("added_loose_norm.root");
c1->SaveAs("added_fr.root");
//c1->SaveAs("added_nt2.root");
//c1->SaveAs("added_nt3_norm.root");
}