#include "drawAsymSysErr.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include <string>
#include <sstream>
#include <fstream>

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

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

double drawBeamAsymmetry(const int maxAsyNum, TH1D* hPERPSub, TH1D* hPARASub, int bin, double norm, double minx, double maxx, int color,  TString figName, TString histName_h)
{
  TString rebin_s = "rebin=" + IntToString(100/bin);
  
  TH1D* hPPSum = new TH1D("hPPSum", histName_h + "_sum", bin, minx, maxx);
  hPPSum->Sumw2();
  hPPSum->Add(hPERPSub, hPARASub, 1, norm); 

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

  cout<<"maximum= "<<maximum<<endl;

  double maximumSum = hPPSum->GetMaximum();
  maximumSum *= 3.5;
  hPPSum->SetMinimum(0);
  hPPSum->SetMaximum(maximumSum);

  // make asymmetry
  TH1I* hAsym = hPERPSub->GetAsymmetry(hPARASub, norm);
  double maximumAsym = hAsym->GetMaximum();
  maximumAsym *= 4.0; 
  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(color);
  hPERPSub->GetYaxis()->SetTitle("Counts/0.314");
  hPERPSub->Fit(fitPERP, "");
  hPERPSub->SetMarkerStyle(20);
  hPERPSub->SetMarkerColor(color);
  hPERPSub->Draw("e");
  //NPERP_btmp = hPERPSub->Integral();
  NPERP_btmp = fitPERP->Integral(minx, maxx);
  //NPERP_btmp = fitPERP->GetNpar();
  TPavesText *paves = new TPavesText(-2.75, 0.8*maximum, -1.4, 0.9*maximum);
  paves->AddText(rebin_s);
  paves->SetFillColor(19);
  paves->Draw();
  
  cc->cd(2);
  TF1* fitPARA = new TF1("psiFitPARA", "[0]*(1.0 + [1]*cos(2*(x + [2])))");
  fitPARA->SetLineColor(color);
  hPARASub->GetYaxis()->SetTitle("Counts/0.314");
  hPARASub->Fit(fitPARA, "");
  hPARASub->Draw("e");
  hPARASub->SetMarkerStyle(20);
  hPARASub->SetMarkerColor(color);
  //NPARA_btmp = fitPARA->GetNpar();
  NPARA_btmp = fitPARA->Integral(minx, maxx);
  //NPARA_btmp = hPARASub->Integral();
        								 
  cc->cd(3);
  TF1* fitAsym = new TF1("asymFit", "[0]*cos(2*(x + [1]))");
  fitAsym->SetLineColor(color);
  hAsym->GetYaxis()->SetTitle("Asymmetry");
  hAsym->Fit(fitAsym, "");
  hAsym->SetMarkerStyle(20);
  hAsym->SetMarkerColor(color);
  hAsym->Draw("e");
  cout<<histName_h<<endl;
  cout<<"PSigma= "<<fitAsym->GetParameter(0)<<" +/- "<<fitAsym->GetParError(0)<<";"<<endl;
  Asymtmp = fitAsym->GetParameter(0);
  Errtmp  = fitAsym->GetParError(0);

  cc->cd(4);
  TF1* fitSum = new TF1("sumFit", "[0]*(1.0 + [1]*cos(2*(x + [2])))");
  fitSum->SetLineColor(color);
  hPPSum->GetYaxis()->SetTitle("Counts/0.314");
  hPPSum->Fit(fitSum, "");
  hPPSum->SetMarkerStyle(20);
  hPPSum->SetMarkerColor(color);
  hPPSum->Draw("e");
   
  cc->Print(figName+".pdf");
  cc->Print(figName+".ps");

  delete hPPSum;
  delete fitPERP;
  delete fitPARA;
  delete fitAsym;
  delete fitSum;
  delete paves;
  delete cc;

  return Asymtmp;
}	

void drawAsymSysErr() {

  //int rebin = 5;
  //int bin   = 100/rebin;
  Double_t normC = 1.0/6.0;
  double norm = 1.10422;  //for ana & systematic err 1.1656
  //double norm = 22022.0/18894.0;  //for ana & systematic err 1.1656
  //double norm = 22112.0/19040.0;  //for systematic err 1.1613
  //double norm = 22081.0/18768.0; //standard  1.1765

  double  minx = -3.14;
  double  maxx = 3.14;

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

  TString headfileName = "./histroot/ana_p2gamma_phip";
  TString histName_head = "Phi_Proton_"; // add any TH1D histogram name
  TString tailfileName = ".root";

  TString figsName_head = "./figPhiP/asym_";

  const int nsysSc=8;  //The number of the source of systematic errors, such as accidentals, coherent beam energy, cut conditions
  const int ncutsRg=3; //The number of the different cut ranges of one source of systematic errors
  const int nPS=2;
  const int maxAsyNum = 3;
  const int nrebin=5;
  //const int nrebin=10;

  TString titleSc[nsysSc]  = {"standard_", "accidentals_", "coherent BE_", "MMsq_", "Delta_phi_", "ME_", "MMXp_", "IM_2gamma_"}; 

  TString titleRg[ncutsRg] = {"cuta_", "cutb_", "cutc_"};

  TString titlePS[nPS] = {"pi0_", "eta_"};

  TString figsName_sec[maxAsyNum] = {"Tot", "Ran", "TSubR"};

  TString title_rebin[nrebin] = {"01_", "02_", "04_", "05_", "10_"};
  int rebin_i[nrebin] = {1, 2, 4, 5, 10};
  //TString title_rebin[nrebin] = {"01_", "02_", "03_", "04_", "05_", "06_", "07_", "08_", "09_", "10_"};

  int color[maxAsyNum] = {1, 4, 2};	

  //for(int rebin = 1; rebin<11; rebin++)
  for(int m = 1; m<nrebin+1; m++)
   {
   int rebin = rebin_i[m-1];
   int bin   = 100/rebin;
   for(int i=0; i<nsysSc; i++)
    {
    for(int j=0; j<ncutsRg; j++)
     {   
     for(int k=0; k<nPS; k++)
      {   
      TString figsName_mid = IntToString(i)+ "_" + IntToString(j)+ "_" + IntToString(k);
      TString perpfileName = headfileName + "_perp" + tailfileName;
      TString parafileName = headfileName + "_para" + tailfileName;

      TString figName[maxAsyNum];
      TString histName = histName_head + figsName_mid;
      TString histName_h;

      TH1D* hPERPlist[maxAsyNum];
      TH1D* hPARAlist[maxAsyNum];


      for(int l=0; l< maxAsyNum; l++)
	 {
         if(i==1&&(j==0||j==1))
            normC = 1.0/8.0;
         else if(i==1&&j==2)
            normC = 1.0/10.0;
         else
            normC = 1.0/6.0;

         histName_h = histName_head + "rebin_" + title_rebin[m-1] + titleSc[i] + titleRg[j] + titlePS[k] + figsName_sec[l];
         hN[m-1][i][j][k][l] = histName_h;
         figName[l]   = figsName_head + "rebin_" + title_rebin[m-1] + titleSc[i] + titleRg[j] + titlePS[k] + figsName_sec[l];

         Errtmp = 0.0;

         TFile *fPERP = TFile::Open(perpfileName);
         if(!fPERP) return;
         TH1D* hPERP = (TH1D*)fPERP->Get(histName);
         hPERP->Sumw2();
         TH1D* hPERPBkg = (TH1D*)fPERP->Get("Acci_" + histName);
         hPERPBkg->Sumw2();
         hPERPBkg->Scale(normC);
         TH1D* hPERPSub = new TH1D("hPERPSub", histName + "_sub", 100, minx, maxx);
         hPERPSub->Sumw2();
         hPERPSub->Add(hPERP, hPERPBkg, 1, -1); 
         hPERPSub->Draw();

         TFile *fPARA = TFile::Open(parafileName);
         if(!fPARA) return;
         TH1D* hPARA = (TH1D*)fPARA->Get(histName);
         hPARA->Sumw2();
         TH1D* hPARABkg = (TH1D*)fPARA->Get("Acci_" + histName);
         hPARABkg->Sumw2();
         hPARABkg->Scale(normC);
         TH1D* hPARASub = new TH1D("hPARASub", histName + "_sub", 100, minx, maxx);
         hPARASub->Sumw2();
         hPARASub->Add(hPARA, hPARABkg, 1, -1);
	 hPERPlist[0]=hPERP;
	 hPERPlist[1]=hPERPBkg;
	 hPERPlist[2]=hPERPSub;
	 hPARAlist[0]=hPARA;
	 hPARAlist[1]=hPARABkg;
	 hPARAlist[2]=hPARASub;

         hPERPlist[l]->Rebin(rebin);
         hPARAlist[l]->Rebin(rebin);
         NPERP_a[m-1][i][j][k][l] = hPERPlist[l]->Integral();
         NPARA_a[m-1][i][j][k][l] = hPARAlist[l]->Integral();
         Norm_a[m-1][i][j][k][l]  = NPERP_a[m-1][i][j][k][l]/NPARA_a[m-1][i][j][k][l];
         Asym[m-1][i][j][k][l] = drawBeamAsymmetry(maxAsyNum, hPERPlist[l], hPARAlist[l], bin, norm, minx, maxx, color[l], figName[l], histName_h);

         Err[m-1][i][j][k][l] = Errtmp;       
         NPERP_b[m-1][i][j][k][l] = NPERP_btmp;
         NPARA_b[m-1][i][j][k][l] = NPARA_btmp;
         Norm_b[m-1][i][j][k][l] = NPERP_b[m-1][i][j][k][l]/NPARA_b[m-1][i][j][k][l];

         delete hPERPSub;
         delete hPARASub;

         fPERP->Close();
         fPARA->Close();
         }
        }
      }
     } 
    }

/*
  for(int m = 1; m<nrebin+1; m++)
   {
   for(int i=0; i<nsysSc; i++)
    {
     for(int j=0; j<ncutsRg; j++)
      {   
       for(int k=0; k<nPS; k++)
        {   
        for(int l=0; l< maxAsyNum; l++)
	 {
         cout<<hN[m-1][i][j][k][l]<<endl;
         cout<<"Asym="<<Asym[m-1][i][j][k][l]<<" +/- "<<Err[m-1][i][j][k][l]<<";"<<endl;
         }
        }
       }
     } 
    }
*/
  cout<<"Beam Asymmetry for pi0 (rebin =5)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {   
    cout<<Asym[3][i][j][0][2]<<" +/- "<<Err[4][i][j][0][2]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

  cout<<"Beam Asymmetry for eta (rebin =5)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {   
    cout<<Asym[3][i][j][1][2]<<" +/- "<<Err[4][i][j][1][2]<<";      ";
    
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

   cout<<"Beam Asymmetry for pi0 with rebin=1, 2, 4, 5, 10 for standard cut condition"<<endl;
   for(int m = 1; m<nrebin+1; m++)
     {
     cout<<Asym[m-1][0][0][0][2]<<" +/- "<<Err[m-1][0][0][0][2]<<";      "<<endl;
     }

   cout<<"Beam Asymmetry for eta with rebin=1, 2, 4, 5, 10 for standard cut condition"<<endl;
   for(int m = 1; m<nrebin+1; m++)
     {
     cout<<Asym[m-1][0][0][1][2]<<" +/- "<<Err[m-1][0][0][1][2]<<";      "<<endl;
     }
/*
  cout<<"The norm from Histograms"<<endl;
  for(int m = 1; m<nrebin+1; m++)
   {
   for(int i=0; i<nsysSc; i++)
    {
     for(int j=0; j<ncutsRg; j++)
      {   
       for(int k=0; k<nPS; k++)
        {   
        for(int l=0; l< maxAsyNum; l++)
	 {
         cout<<hN[m-1][i][j][k][l]<<endl;
         cout<<"NPERP_a="<<NPERP_a[m-1][i][j][k][l]<<";   NPERP_b="<<NPERP_b[m-1][i][j][k][l]<<";"<<endl;
         cout<<"NPARA_a="<<NPARA_a[m-1][i][j][k][l]<<";   NPARA_b="<<NPARA_b[m-1][i][j][k][l]<<";"<<endl;
         cout<<"norm_a="<<Norm_a[m-1][i][j][k][l]<<";   norm_b="<<Norm_b[m-1][i][j][k][l]<<";"<<endl;
         }
        }
       }
     } 
    }
  */

  return;
}