#include "TFile.h"
#include "TTree.h"
#include "TH2F.h"

Int_t PlotKTUY() {
  TFile *f1 = new TFile("KTUY05_m246.root","READ");
  TTree *tree1 = (TTree*)f1->Get("KTUY05_m246");
  TH2F  *hist = new TH2F("hist", "KTUY05 Mass Excess (MeV);Nutrom number N; Atomic number Z;",201,-0.5,200.5,131,-0.5,130.5);
  TH2F  *hist2= new TH2F("hist2","",201,-0.5,200.5,131,-0.5,130.5);
  tree1->Draw("ZZ:NN>>hist","Mcal*(bBound ==1)");
  tree1->Draw("ZZ:NN>>hist2","bBound == 1");
  
  BGSub(hist,hist2);
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->CenterTitle();
  hist->Draw("COLZ");
  
  //  delete hist;
  //  delete f;
  return 0;
}

Int_t BGSub(TH2F *hist,TH2F *hist2) {
  Double_t zmin;
  TH2F* hist3 = hist->Clone("hist3");
  Fill_1(hist3);
  hist3->Add(hist3,hist2,1.0,-1.0);
  zmin = hist->GetMinimum();
  hist->Add(hist,hist3,1.0,zmin*1.1);
  hist->SetMinimum(zmin);
  return 0;
}

Int_t Fill_1(TH2F *hist) {
  Int_t xbins,ybins;
  Double_t xcnt,ycnt;

  xbins = hist->GetNbinsX();
  ybins = hist->GetNbinsY();

  hist->Reset();
  for (Int_t i=1; i <= xbins; i++) {
    for (Int_t j=1; j <= ybins; j++) {
      xcnt = hist->GetXaxis()->GetBinCenter(i);
      ycnt = hist->GetYaxis()->GetBinCenter(j);
      hist->Fill(xcnt,ycnt);
    }
  }
  return 0;
}