#include "drawAsym.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"
#include "TAxis.h"

double integ(TH1D* histdata, double xmin, double xmax)
{
  TAxis *axis = histdata->GetXaxis();
  int bmin = axis->FindBin(xmin);
  int bmax = axis->FindBin(xmax);
  double para = histdata->Integral(bmin, bmax);
  return para;
}

void drawTH1D2Compare(TH1D* hist1, TH1D* hist2, TString figName)
{
  Int_t w = 960;
  Int_t h = 600; 
  double minimum = 0.1;
  double maximum = 50*hist1->GetMaximum();
  TString title_c = "Im_2g_omega";
  TCanvas *c = new TCanvas("c", title_c, w, h);
  c->Divide(1,1);
  hist1->SetMinimum(minimum);
  gPad->SetLogy();
  hist1->SetMaximum(maximum);
  hist1->SetLineColor(kRed);
  hist1->SetLineStyle(1);
  hist2->SetLineColor(kBlue);
  hist2->SetLineStyle(2);
  hist1->Draw(); 
  hist2->Draw("SAME");

  TLegend *leg = new TLegend(0.42,0.72, 0.72, 0.86);
  leg->SetHeader("The omega background");
  leg->SetTextSize(0.04);
  leg->AddEntry(hist1, "data", "f");
  leg->AddEntry(hist2, "#omega mc", "f");
  leg->Draw();
  c->Print(figName);
  return;
}

void drawomega()
{
  TString fileNamedata = "./histroot/ana_p2gamma_Cuts_NewData.root";
  TFile *fdata = TFile::Open(fileNamedata);
  TString fileNameinmc = "./histroot/hist_p2gamma_inc_010000.root";
  TFile *finmc = TFile::Open(fileNameinmc);
  TString fileNameomega = "./histroot/hist_p2gamma_omega_010000.root";
  TFile *fomega = TFile::Open(fileNameomega);
  
  TString figsDir = "./figsomega/"; 
    cout<<"root file is: "<<fileNamedata.Data()<<endl;
    cout<<"root file is: "<<fileNameinmc.Data()<<endl;
    cout<<"root file is: "<<fileNameomega.Data()<<endl;

  double xmin = 0.660; 
  double xmax = 0.900;
  double eta_min = 0.482;
  double eta_max = 0.610;

  TString histName = "S_IM_2gamma_1";
  TString figName  = figsDir + "Im_2g_omega.pdf";
  TH1D* histdata   = (TH1D*)fdata->Get(histName);
  TH1D* histomega  = (TH1D*)fomega->Get(histName);
  double paradata  = integ(histdata, xmin, xmax);
  double paraomega = integ(histomega, xmin, xmax);
  double paraomega_eta_bf = integ(histomega, eta_min, eta_max);
  cout<<"paraomega_eta_bf = "<<paraomega_eta_bf<<endl;
  double norm = paradata/paraomega;
  histomega->Scale(norm);
  double paraomega_eta = integ(histomega, eta_min, eta_max); 
  double paradata_eta = integ(histdata, eta_min, eta_max); 
  double fraction = paraomega_eta/paradata_eta;
  cout<<"paraomega_eta_af = "<<paraomega_eta<<endl;
  cout<<"fraction = "<<fraction<<endl;
  drawTH1D2Compare(histdata, histomega, figName);
  double PStot=0.3519;
  double ff = 0.0038;
  double PSbkg = -0.1736;
  //double PStot = Asym[0][1][2];
  //double PSbkg = Asymomega; 
  //double ff = fraction;
  double PSsig = (PStot-ff*PSbkg)/(1-ff);
  cout<<Asym[0][1][2]<<endl;
  cout<<"PSsig(eta not corrected)="<<PStot<<endl;
  cout<<"PSsig(eta corrected)="<<PSsig<<endl;
}