#include "drawAsym.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();
}

void drawBeamAsymmetry(const int maxAsyNum, TH1D* hPERPSub, TH1D* hPARASub, int i, int j, int k, int bin, double norm, double minx, double maxx, TString figName)
{
  int color[maxAsyNum] = {1, 4, 2};	
  const int maxhistNum = 2;
  TString histName[maxhistNum] = {"Phi_Proton_pi0", "Phi_Proton_eta"}; // add any TH1D histogram name
  TH1D* hPPSum = new TH1D(Form("hPERPSum_%d_%d_%d", i, j, k), histName[j] + "_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.0;
  hPPSum->SetMinimum(0);
  hPPSum->SetMaximum(maximumSum);

  // make asymmetry
  TH1I* hAsym = hPERPSub->GetAsymmetry(hPARASub, norm);
  double maximumAsym = hAsym->GetMaximum();
  maximumAsym *= 3.0; 
  hAsym->SetMaximum(maximumAsym);
  
  TCanvas *cc = new TCanvas(Form("cc_%d_%d_%d", i, j, k), "cc", 900, 760);
  cc->Divide(2,2);
  cc->cd(1);
  TF1* fitPERP = new TF1(Form("psiFitPERP_%d", k), "[0]*(1.0 + [1]*cos(2*(x + [2])))");
  fitPERP->SetLineColor(color[k]);
  hPERPSub->GetYaxis()->SetTitle("Counts/0.314");
  hPERPSub->Fit(fitPERP, "");
  hPERPSub->SetMarkerStyle(20);
  hPERPSub->SetMarkerColor(color[k]);
  hPERPSub->Draw("e");
  
  cc->cd(2);
  TF1* fitPARA = new TF1(Form("psiFitPARA_%d", k), "[0]*(1.0 + [1]*cos(2*(x + [2])))");
  fitPARA->SetLineColor(color[k]);
  hPARASub->GetYaxis()->SetTitle("Counts/0.314");
  hPARASub->Fit(fitPARA, "");
  hPARASub->Draw("e");
  hPARASub->SetMarkerStyle(20);
  hPARASub->SetMarkerColor(color[k]);
        								 
  cc->cd(3);
  TF1* fitAsym = new TF1(Form("asymFit_%d", k), "[0]*cos(2*(x + [1]))");
  fitAsym->SetLineColor(color[k]);
  hAsym->GetYaxis()->SetTitle("Asymmetry");
  hAsym->Fit(fitAsym, "");
  hAsym->SetMarkerStyle(20);
  hAsym->SetMarkerColor(color[k]);
  hAsym->Draw("e");
  if(i==0&&j==0&&k==2)
     cout<<"PSigma(pi0)= "<<fitAsym->GetParameter(0)<<" +/- "<<fitAsym->GetParError(0)<<";"<<endl;
  if(i==0&&j==1&&k==2)
     cout<<"PSigma(eta)= "<<fitAsym->GetParameter(0)<<" +/- "<<fitAsym->GetParError(0)<<";"<<endl;
  Asym[i][j][k] = fitAsym->GetParameter(0);

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

void drawAsym() {

  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;

  double  minx = -3.14;
  double  maxx = 3.14;

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

  const int maxdatafileNum = 1;
  const int maxhistNum = 2;
  TString headfileName = "./histroot/ana_p2gamma_phip";
  TString midfileName[maxdatafileNum] = {""}; // add any TH1D histogram name
  //TString midfileName[maxdatafileNum] = {"", "_1049", "_107", "_1087", "_1083"}; // add any TH1D histogram name
  TString histName[maxhistNum] = {"Phi_Proton_pi0", "Phi_Proton_eta"}; // add any TH1D histogram name
  TString figsubName[maxhistNum] = {"pi0_", "eta_"}; // add any TH1D histogram name
  TString tailfileName = ".root";

  TString figsName_head = "./figPhiP/asym_";
  const int maxAsyNum = 3;
  TString figsName_sec[maxAsyNum] = {"Tot", "Ran", "TSubR"};

  for(int i=0; i<maxdatafileNum; i++)
    {
    for(int j=0; j<maxhistNum; j++)
      {   
      string figsName_mid = IntToString(i)+ "_" + IntToString(j);
      TString perpfileName = headfileName + "_perp" + midfileName[i] + tailfileName;
      TString parafileName = headfileName + "_para" + midfileName[i] + tailfileName;

      TString figName[maxAsyNum];

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

      for(int k=0; k< maxAsyNum; k++)
	 {
         TFile *fPERP = TFile::Open(perpfileName);
         if(!fPERP) return;
         TH1D* hPERP = (TH1D*)fPERP->Get(histName[j]);
         hPERP->Sumw2();
         TH1D* hPERPBkg = (TH1D*)fPERP->Get("Acci_" + histName[j]);
         hPERPBkg->Sumw2();
         hPERPBkg->Scale(normC);
         TH1D* hPERPSub = new TH1D(Form("hPERPSub_%d_%d_%d", i, j, k), histName[j] + "_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[j]);
         hPARA->Sumw2();
         TH1D* hPARABkg = (TH1D*)fPARA->Get("Acci_" + histName[j]);
         hPARABkg->Sumw2();
         hPARABkg->Scale(normC);
         TH1D* hPARASub = new TH1D(Form("hPARASub_%d_%d_%d", i, j, k), histName[j] + "_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[k]->Rebin(rebin);
         hPARAlist[k]->Rebin(rebin);
         figName[k]   = figsName_head + figsubName[j] + figsName_sec[k] + midfileName[i];
         drawBeamAsymmetry(maxAsyNum, hPERPlist[k], hPARAlist[k], i, j, k, bin, norm, minx, maxx, figName[k]);
         }
       } 
     }
  cout<<"Asym of pi0="<<Asym[0][0][2]<<endl;
  cout<<"Asym of eta="<<Asym[0][1][2]<<endl;

  return;
}