#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
#include <string>
#include <sstream>
#include "TGraphErrors.h"

string DoubleToString(double a)
{
  ostringstream temp;
  temp<<a;
  return temp.str();
}

string IntToString(int a)
{
  ostringstream temp;
  temp<<a;
  return temp.str();
}

void drawbinnedPSigma(double x[], double y[], double ey[], const int n, TString figName)
{
   gStyle->SetOptStat(0);

   TCanvas *c1 = new TCanvas("c1", "The centroid with error bars", 100, 5, 1000, 600);
   c1->Divide(1,1); 

   double ex[n] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
   c1->cd(1);
   TH1F* h = new TH1F("h", "h", 100, 0., 2.);
   h->Draw();
   TGraphErrors *gr = new TGraphErrors(n, x, y, ex, ey);
   //gr->Draw("AP*");
   h->SetMinimum(0.10); 
   h->SetMaximum(0.50);
   gr->SetMarkerStyle(20);
   gr->SetMarkerSize(1);
   gr->SetMarkerColor(kBlack);
   gr->SetLineColor(kBlack);
   h->GetXaxis()->SetTitle("<-t> (GeV^{2})"); 
   h->GetYaxis()->SetTitle("P #Sigma"); 
   h->GetXaxis()->SetRangeUser(0.0,2.0);
   h->GetYaxis()->SetTitleOffset(1.2);
   h->SetTitle("The P#Sigma of #eta^{0} photoproduction for total");
   //c1->Update();
   gr->Draw("same P");
   //c1->Print("figPhiP/PSigma.pdf");
   c1->Print(figName);
   return;

} 
 
double Meanmt(TString datafileName, TString Hist2DName, double minx, double maxx, TString figsName, int i)
{
  TFile *fdata = TFile::Open(datafileName);
  TH2D* h = (TH2D*)fdata->Get(Hist2DName);

  gStyle->SetOptStat("nemr");

  //Double_t ylow = 0.0;
  //Double_t yup = 0.2;

  Double_t ylow = -3.14;
  Double_t yup = 3.14;

  TCanvas *cc_ave = new TCanvas("cc_ave", "cc_ave", 600, 400);
  cc_ave->Divide(1,1);
  double maximum=0.0;
  int binMiny = h->GetYaxis()->FindBin(ylow);
  int binMaxy = h->GetYaxis()->FindBin(yup);
  TH1D *h1 = (TH1D*)h->ProjectionX(Form("E2gamma_%d", i),binMiny,binMaxy);
  maximum = h1->GetMaximum();
  maximum *= 3.6;
  h1->SetMaximum(maximum);

  int binMinx = h->GetXaxis()->FindBin(minx);
  int binMaxx = h->GetXaxis()->FindBin(maxx);
  h1->GetXaxis()->SetRange(binMinx,binMaxx);
  cout<<"min_x="<<h1->GetMean(1)<<endl;
  cc_ave->cd(1);
  h1->Draw("");
  cc_ave->Print(figsName);
  return h1->GetMean(1);
}




void drawetaAsysBkgTSubRbin(double lowCoh = 8.4, double highCoh = 9.0) {

  int rebin = 5;

  double  minx = -3.14;
  double  maxx = 3.14;

  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1111);

  TFile *fPERP = TFile::Open("./histroot/ana_p2gamma_phip_perp.root");
  if(!fPERP) return;
  TH2D* h2PERP = (TH2D*)fPERP->Get("Phi_Proton_vs_t_5");
  TH2D* h2PERPBkg = (TH2D*)fPERP->Get("Acci_Phi_Proton_vs_t_5");

  TFile *fPARA = TFile::Open("./histroot/ana_p2gamma_phip_para.root");
  if(!fPARA) return;
  TH2D* h2PARA = (TH2D*)fPARA->Get("Phi_Proton_vs_t_5");
  TH2D* h2PARABkg = (TH2D*)fPARA->Get("Acci_Phi_Proton_vs_t_5");

  TString datafileName = "./histroot/ana_p2gamma_Cuts_NewData.root";
  TString Hist2DName = "Phi_Proton_vs_t_5";

  const int nt = 7;
  double Array_t[nt] = {0.00, 0.14, 0.20, 0.30, 0.68, 1.0, 2.0};
  TString Array_t_s[nt];
  for(int i = 0; i<nt; i++)
     {
     Array_t_s[i] = DoubleToString(Array_t[i]);
     }

  double x[nt-1]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  double y[nt-1]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  double ey[nt-1]={0.0, 0.0, 0.0, 0.0, 0.0, 0.0};


  cout<<"Array_t[0]="<<Array_t[0]<<endl;
  cout<<"Array_t[1]="<<Array_t[1]<<endl;
  cout<<"Array_t[2]="<<Array_t[2]<<endl;

  cout<<"Array_t_s[0]="<<Array_t_s[0]<<endl;
  cout<<"Array_t_s[1]="<<Array_t_s[1]<<endl;
  cout<<"Array_t_s[2]="<<Array_t_s[2]<<endl;

  string figsName_head = "./figPhiP/asym_tbin_etaTSubR_";
  string figsName_tail = ".pdf";

  //for(int i=1; i<(nt); i++)
  for(int i=1; i<7; i++)
     {
     string figsName_mid = IntToString(i);
     TString figsName = figsName_head + figsName_mid + figsName_tail;
     TString figsName1 = figsName_head + figsName_mid + "test" +  figsName_tail;
     TString figsmtName = "./figPhiP/mt" + figsName_mid + figsName_tail;

     int j=i-1;
     cout<<"j="<<j<<endl;
     cout<<"i="<<i<<endl;
     cout<<"Array_t[j]="<<Array_t[j]<<endl;
     cout<<"Array_t[i]="<<Array_t[i]<<endl;

     x[j] = Meanmt(datafileName, Hist2DName, Array_t[j], Array_t[i], figsmtName, i);

     int binMinPERP = h2PERP->GetXaxis()->FindBin(Array_t[j]);
     int binMaxPERP = h2PERP->GetXaxis()->FindBin(Array_t[i]);
     TH1D *hPERP = (TH1D*)h2PERP->ProjectionY(Form("phi_p_perp_%d", i),binMinPERP,binMaxPERP);
     hPERP->Sumw2();
     hPERP->Rebin(rebin);
     int binMinPERPBkg = h2PERPBkg->GetXaxis()->FindBin(Array_t[j]);
     int binMaxPERPBkg = h2PERPBkg->GetXaxis()->FindBin(Array_t[i]);
     TH1D *hPERPBkg = (TH1D*)h2PERPBkg->ProjectionY(Form("phi_p_perp_bkg_%d", i),binMinPERPBkg,binMaxPERPBkg);
     hPERPBkg->Sumw2();
     hPERPBkg->Rebin(rebin);
     Double_t normC = 1.0/6.0;
     hPERPBkg->Scale(normC);
     TH1D* hPERPSub = new TH1D("hPERPSub", "Phi_Proton_eta_cutBE_sub", 20, minx, maxx);
     hPERPSub->Sumw2();
     hPERPSub->Add(hPERP, hPERPBkg, 1, -1); 
     hPERPSub->Draw();


     int binMinPARA = h2PARA->GetXaxis()->FindBin(Array_t[j]);
     int binMaxPARA = h2PARA->GetXaxis()->FindBin(Array_t[i]);
     TH1D *hPARA = (TH1D*)h2PARA->ProjectionY(Form("phi_p_para_%d", i),binMinPARA,binMaxPARA);
     hPARA->Sumw2();
     hPARA->Rebin(rebin);
     int binMinPARABkg = h2PARABkg->GetXaxis()->FindBin(Array_t[j]);
     int binMaxPARABkg = h2PARABkg->GetXaxis()->FindBin(Array_t[i]);
     TH1D *hPARABkg = (TH1D*)h2PARABkg->ProjectionY(Form("phi_p_para_bkg_%d", i),binMinPARABkg,binMaxPARABkg);
     hPARABkg->Sumw2();
     hPARABkg->Rebin(rebin);
     hPARABkg->Scale(normC);
     TH1D* hPARASub = new TH1D("hPARASub", "Phi_Proton_eta_cutBE_sub", 20, minx, maxx);
     hPARASub->Sumw2();
     hPARASub->Add(hPARA, hPARABkg, 1, -1); 

     //double norm = hPERPSub->Integral()/hPARASub->Integral();
     double norm = 22081.0/18768.0;
     //hPARASub->Scale(norm);

     // set max-min and get normalization
     double maximum = hPERPSub->GetMaximum();
     if(hPARASub->GetMaximum() > maximum)
       maximum = hPARA->GetMaximum();
     maximum *= 2.6;
     hPERPSub->SetMinimum(0);
     hPARASub->SetMinimum(0);
     hPERPSub->SetMaximum(maximum);
     hPARASub->SetMaximum(maximum);

     TH1D* hPPSum = new TH1D("hPERPSum", "Phi_Proton_eta_cutBE_sum", 20, minx, maxx);
     hPPSum->Sumw2();
     hPPSum->Add(hPERPSub, hPARASub, 1, norm); 
     
     double maximumSum = hPPSum->GetMaximum();
     maximumSum *= 2.6;
     hPPSum->SetMinimum(0);
     hPPSum->SetMaximum(maximumSum);

     // make asymmetry
     TH1I* hAsym = hPERPSub->GetAsymmetry(hPARASub, norm);
     double maximumAsym = hAsym->GetMaximum();
     maximumAsym *= 2.6; 
     hAsym->SetMaximum(maximumAsym);
     
     TCanvas *cc = new TCanvas("cc", "cc", 900, 760);
     cc->Divide(2,2);
     cc->cd(1);
     TF1* fitPERP = new TF1("psiFitPERP","[0]*(1.0 + [1]*cos(2*(x + [2])))");
     fitPERP->SetLineColor(kRed);
     hPERPSub->Fit(fitPERP, "");
     hPERPSub->SetMarkerStyle(20);
     hPERPSub->SetMarkerColor(kRed);
     hPERPSub->Draw("e");
     
     cc->cd(2);
     TF1* fitPARA = new TF1("psiFitPARA","[0]*(1.0 + [1]*cos(2*(x + [2])))");
     fitPARA->SetLineColor(kRed);
     hPARASub->Fit(fitPARA, "");
     hPARASub->Draw("e");
     hPARASub->SetMarkerStyle(20);
     hPARASub->SetMarkerColor(kRed);
           								 
     cc->cd(3);
     TF1* fitAsym = new TF1("asymFit","[0]*cos(2*(x + [1]))");
     fitAsym->SetLineColor(kRed);
     hAsym->Fit(fitAsym, "");
     hAsym->SetMarkerStyle(20);
     hAsym->SetMarkerColor(kRed);
     hAsym->Draw("e");
     y[j]  = fitAsym->GetParameter(0);
     ey[j] = fitAsym->GetParError(0);

     y[j] = (y[j]+0.0038*0.1736)/(1-0.0038);

     cc->cd(4);
     TF1* fitSum = new TF1("sumFit","[0]*(1.0 + [1]*cos(2*(x + [2])))");
     fitSum->SetLineColor(kRed);
     hPPSum->Fit(fitSum, "");
     hPPSum->SetMarkerStyle(20);
     hPPSum->SetMarkerColor(kRed);
     hPPSum->Draw("e");

  
     cc->Print(figsName);

     cout<<"x[j]="<<x[j]<<"; y[j]="<<y[j]<<"; ey[j]="<<ey[j]<<"; nt="<<nt<<endl;
 }

  drawbinnedPSigma(x, y, ey, nt, "figPhiP/etaPSigma.pdf");
 
  return;
}