#include "drawAsymSysErrOme.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_bin(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");
  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);
        								 
  cc->cd(3);
  TF1* fitAsym = new TF1("asymFit", "[0]*cos(2*(x + [1]))");
  fitAsym->FixParameter(1, -0.07774);
  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->FixParameter(2, -0.07774);
  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 drawbinnedPSigma3_sys(double xx[][nsysSc][3][6], double yy[][nsysSc][3][6], double eyy[][nsysSc][3][6], const int n, const int nrebin, const int nsysSc, const int ncutsRg, const int nPS, TString figName)
{
   gStyle->SetOptStat(0);

   double x[5][nsysSc][3][3], y[5][nsysSc][3][3], ey[5][nsysSc][3][3];

   double exx[n] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
   double ex[3] = {0.0, 0.0, 0.0};
  for(int m = 1; m<nrebin+1; m++)
   {
   for(int i=0; i<nsysSc; i++)
    {
    for(int j=0; j<ncutsRg; j++)
     {   
      for(int ii=0; ii<3; ii++)
         {
         x[m-1][i][j][ii]=xx[m-1][i][j][2*ii+1];
         y[m-1][i][j][ii]=yy[m-1][i][j][2*ii+1];
         ex[ii]=exx[2*ii+1];
         ey[m-1][i][j][ii]=eyy[m-1][i][j][2*ii+1];
         }
       }
      }
     }
      

   double x_0[n], x_1[n], x_2[n], Asy_0[n], Asy_1[n], Asy_2[n], Err_0[n], Err_1[n], Err_2[n];
   int rebin_i[nrebin] = {1, 2, 4, 5, 10};
   int bin;
   double x_bins_0[n], x_bins_1[n], x_bins_2[n], x_bins_3[n], x_bins_4[n],  Asy_bins_0[n], Asy_bins_1[n], Asy_bins_2[n], Asy_bins_3[n], Asy_bins_4[n], Err_bins_0[n], Err_bins_1[n], Err_bins_2[n], Err_bins_3[n], Err_bins_4[n];

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

    if(1) 
      {

       TString fig_bins_index = "_bins_m2g";

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

       c1->cd(1);
       TH1F* h = new TH1F("h", "h", 100, 0., 2.0);
       h->Draw();
       h->SetMinimum(-0.60); 
       h->SetMaximum(0.80);
       h->GetXaxis()->SetTitle("m_{2#gamma} (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 photoproduction for TSubR");

       TString rebin_s = "different bins standard cut";
       TString note_color = "black for 100, blue for 50, red for 25, green for 20, yellow for 10";
       TPavesText *paves = new TPavesText(0.4, 0.7, 1.2, 0.9);
       paves->AddText(rebin_s);
       paves->AddText(note_color);
       paves->SetFillColor(19);
       paves->Draw();
       
      for(int l=0; l<3; l++)
       {
  cout<<__LINE__<<endl;
       x_bins_0[l]   = x[0][0][0][l];
       Asy_bins_0[l] = y[0][0][0][l];
       Err_bins_0[l] = ey[0][0][0][l]; 
       x_bins_1[l]   = x[1][0][0][l];
       Asy_bins_1[l] = y[1][0][0][l];
       Err_bins_1[l] = ey[1][0][0][l]; 
       x_bins_2[l]   = x[2][0][0][l];
       Asy_bins_2[l] = y[2][0][0][l];
       Err_bins_2[l] = ey[2][0][0][l]; 
       x_bins_3[l]   = x[3][0][0][l];
       Asy_bins_3[l] = y[3][0][0][l];
       Err_bins_3[l] = ey[3][0][0][l]; 
       x_bins_4[l]   = x[4][0][0][l];
       Asy_bins_4[l] = y[4][0][0][l];
       Err_bins_4[l] = ey[4][0][0][l]; 
        }   
       TGraphErrors *gr0 = new TGraphErrors(n, x_bins_0, Asy_bins_0, ex, Err_bins_0);
       TGraphErrors *gr1 = new TGraphErrors(n, x_bins_1, Asy_bins_1, ex, Err_bins_1);
       TGraphErrors *gr2 = new TGraphErrors(n, x_bins_2, Asy_bins_2, ex, Err_bins_2);
       TGraphErrors *gr3 = new TGraphErrors(n, x_bins_3, Asy_bins_3, ex, Err_bins_3);
       TGraphErrors *gr4 = new TGraphErrors(n, x_bins_4, Asy_bins_4, ex, Err_bins_4);
       gr0->SetMarkerStyle(20);
       gr0->SetMarkerSize(1);
       gr0->SetMarkerColor(kBlack);
       gr0->SetLineColor(kBlack);
       gr0->Draw("same P");
       gr1->SetMarkerStyle(21);
       gr1->SetMarkerSize(1);
       gr1->SetMarkerColor(kBlue);
       gr1->SetLineColor(kBlue);
       gr1->Draw("same P");
       gr2->SetMarkerStyle(22);
       gr2->SetMarkerSize(1);
       gr2->SetMarkerColor(kRed);
       gr2->SetLineColor(kRed);
       gr2->Draw("same P");
       gr3->SetMarkerStyle(23);
       gr3->SetMarkerSize(1);
       gr3->SetMarkerColor(kGreen);
       gr3->SetLineColor(kGreen);
       gr3->Draw("same P");
       gr4->SetMarkerStyle(24);
       gr4->SetMarkerSize(1);
       gr4->SetMarkerColor(kYellow);
       gr4->SetLineColor(kYellow);
       gr4->Draw("same P");

       c1->Print(figName + fig_bins_index + ".ps");
       c1->Print(figName + fig_bins_index + ".pdf");
  
       delete gr0;
       delete gr1;
       delete gr2;
       delete gr3;
       delete gr4;
       delete h;
       delete c1;
       delete paves;

      } 

  for(int m = 1; m<nrebin+1; m++)
   {
       bin = 100/rebin_i[m-1];
   for(int i=0; i<nsysSc; i++)
    {

       TString fig_index = "_" + IntToString(m-1) + "_" + IntToString(i);
       TString bin_s = IntToString(bin);

       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.0);
       h->Draw();
       h->SetMinimum(-0.60); 
       h->SetMaximum(0.80);
       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 #pi^{0} photoproduction for total");

       TString rebin_s = bin_s + "bins" + titleSc[i];
       TString note_color = "black for cuta, blue for cutb, red for cutc";
       TPavesText *paves = new TPavesText(0.4, 0.7, 1.2, 0.9);
       paves->AddText(rebin_s);
       paves->AddText(note_color);
       paves->SetFillColor(19);
       paves->Draw();
       
      for(int l=0; l<3; l++)
       {
      cout<<y[m-1][i][0][l]<<",   ";
      cout<<ey[m-1][i][0][l]<<";   ";

      if(l==2)  cout<<endl;
  cout<<__LINE__<<endl;
       x_0[l]   = x[m-1][i][0][l];
       Asy_0[l] = y[m-1][i][0][l];
       Err_0[l] = ey[m-1][i][0][l]; 
       x_1[l]   = x[m-1][i][1][l];
       Asy_1[l] = y[m-1][i][1][l];
       Err_1[l] = ey[m-1][i][1][l]; 
       x_2[l]   = x[m-1][i][2][l];
       Asy_2[l] = y[m-1][i][2][l];
       Err_2[l] = ey[m-1][i][2][l]; 
        }   
       TGraphErrors *gr0 = new TGraphErrors(n, x_0, Asy_0, ex, Err_0);
       TGraphErrors *gr1 = new TGraphErrors(n, x_1, Asy_1, ex, Err_1);
       TGraphErrors *gr2 = new TGraphErrors(n, x_2, Asy_2, ex, Err_2);
       gr0->SetMarkerStyle(20);
       gr0->SetMarkerSize(1);
       gr0->SetMarkerColor(kBlack);
       gr0->SetLineColor(kBlack);
       gr0->Draw("same P");
       gr1->SetMarkerStyle(21);
       gr1->SetMarkerSize(1);
       gr1->SetMarkerColor(kBlue);
       gr1->SetLineColor(kBlue);
       gr1->Draw("same P");
       gr2->SetMarkerStyle(22);
       gr2->SetMarkerSize(1);
       gr2->SetMarkerColor(kRed);
       gr2->SetLineColor(kRed);
       gr2->Draw("same P");

  cout<<__LINE__<<endl;
       c1->Print(figName + fig_index + ".ps");
       c1->Print(figName + fig_index + ".pdf");
  cout<<__LINE__<<endl;
  
       delete gr0;
       delete gr1;
       delete gr2;
       delete h;
       delete c1;
       delete paves;

       
      }
   }


   
   return;

} 

/* 
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.0);
   h->Draw();
   TGraphErrors *gr = new TGraphErrors(n, x, y, ex, ey);
   //gr->Draw("AP*");
   h->SetMinimum(0.00); 
   h->SetMaximum(1.00);
   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 #pi^{0} photoproduction for total");
   //c1->Update();
   gr->Draw("same P");
   //c1->Print("figPhiP/PSigma.pdf");
   c1->Print(figName);
  
   delete c1;
   delete h;
   delete gr;
   
   return;

} 
*/


double Meanmt(TString datafileName, TString Hist2DName, double minx, double maxx, int i)
{
  cout<<"beginning of Meanmt"<<endl;
  TFile *fdata = TFile::Open(datafileName);
  TH2D* h = (TH2D*)fdata->Get(Hist2DName);

  cout<<datafileName<<";   "<<Hist2DName<<endl;

  gStyle->SetOptStat("nemr");

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

  Double_t ylow = -3.14;
  Double_t yup = 3.14;

  cout<<__LINE__<<endl;

  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);
  cout<<__LINE__<<endl;

  TH1D *h1 = (TH1D*)h->ProjectionX(Form("E2gamma_%d", i),binMiny,binMaxy);
  cout<<__LINE__<<endl;

  maximum = h1->GetMaximum();
  maximum *= 3.6;
  h1->SetMaximum(maximum);

  cout<<__LINE__<<endl;

  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("");
  double meantmp = h1->GetMean(1);
  //cc_ave->Print(figsName);
  delete cc_ave;
  delete fdata;
  return meantmp;
}




void drawomegaAsymSysErr(double norm) 
{
  
  cout<<"beginning of drawomegaAsymSysErr"<<endl;

  double  minx = -3.14;
  double  maxx = 3.14;
  Double_t normC = 1.0/6.0;

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

  TString hist2DName_head = "Phi_Proton_vs_IM2g_cutBE_sys_";
  TString perpfileName = "./histroot/ana_p2gamma_phip_perp.root";
  TString parafileName = "./histroot/ana_p2gamma_phip_para.root";
   
  TString figsName_head = "./figPhiP/m2gbin/asym_m2gbin_TSubR_";

  cout<<__LINE__<<endl;
  //const int nsysSc=9;  //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 maxAsyNum = 3;
  const int nrebin=5;
  const int nPS = 2;
  const int nt = 7;

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

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

  TString title_rebin[nrebin] = {"01_", "02_", "04_", "05_", "10_"};
  int rebin_i[nrebin] = {1, 2, 4, 5, 10};
  TString titletbin[nt-1] = {"01_", "02_", "03_", "04_", "05_", "06_"};

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

  cout<<"before loop"<<endl;

  //for(int m = 1; m<2; m++)
  for(int m = 1; m<nrebin+1; m++)
   {

   //if(m>2)  continue;

   int rebin = rebin_i[m-1];
   int bin   = 100/rebin;
   for(int i=0; i<nsysSc; i++)
    {

    //if(i>2)  continue;

    for(int j=0; j<ncutsRg; j++)
     {   
         
         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;

      TString histName_mid = IntToString(i)+ "_" + IntToString(j);

      TString figName;
      TString hist2DName = hist2DName_head + histName_mid;
      TString histName_h;


         Errtmp = 0.0;

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

  cout<<"before Array"<<endl;

  const int nt = 7;
  double Array_t[nt] = {0.0, 0.114, 0.156, 0.482, 0.610, 0.660, 0.900};

  TString Array_t_s[nt];
  for(int ii = 0; ii<nt; ii++)
     {
     Array_t_s[ii] = DoubleToString(Array_t[ii]);
     }

  double x[nrebin][nsysSc][ncutsRg][nt-1];
  double y[nt-1];
  double ey[nt-1];


  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;

  
  for(int ii=1; ii<nt; ii++)
     {
     histName_h = hist2DName_head + "rebin_" + title_rebin[m-1] + titleSc[i] + titleRg[j] +titletbin[ii];

     cout<<perpfileName<<endl;
     cout<<parafileName<<endl;
     cout<<hist2DName<<endl;
     cout<<"Acci_" + hist2DName<<endl;

     TFile *fPERP = TFile::Open(perpfileName);
     if(!fPERP) return;
     TH2D* h2PERP = (TH2D*)fPERP->Get(hist2DName);
     TH2D* h2PERPBkg = (TH2D*)fPERP->Get("Acci_" + hist2DName);

     TFile *fPARA = TFile::Open(parafileName);
     if(!fPARA) return;
     TH2D* h2PARA = (TH2D*)fPARA->Get(hist2DName);
     TH2D* h2PARABkg = (TH2D*)fPARA->Get("Acci_" + hist2DName);

     TString figsName_mid = IntToString(ii);
     figName   = figsName_head + "rebin_" + title_rebin[m-1] + titleSc[i] + titleRg[j] + figsName_mid;

     int jj=ii-1;
     cout<<"jj="<<jj<<endl;
     cout<<"ii="<<ii<<endl;
     cout<<"Array_t[jj]="<<Array_t[jj]<<endl;
     cout<<"Array_t[ii]="<<Array_t[ii]<<endl;

     x[m-1][i][j][jj] = Meanmt(datafileName, hist2DName, Array_t[jj], Array_t[ii], ii);
  cout<<__LINE__<<endl;

     int binMinPERP = h2PERP->GetXaxis()->FindBin(Array_t[jj]);
     int binMaxPERP = h2PERP->GetXaxis()->FindBin(Array_t[ii]);
     TH1D *hPERP = (TH1D*)h2PERP->ProjectionY(Form("phi_p_perp_%d", ii),binMinPERP,binMaxPERP);
     hPERP->Sumw2();
     //hPERP->Rebin(rebin);
     int binMinPERPBkg = h2PERPBkg->GetXaxis()->FindBin(Array_t[jj]);
     int binMaxPERPBkg = h2PERPBkg->GetXaxis()->FindBin(Array_t[ii]);
     TH1D *hPERPBkg = (TH1D*)h2PERPBkg->ProjectionY(Form("phi_p_perp_bkg_%d", ii),binMinPERPBkg,binMaxPERPBkg);
     hPERPBkg->Sumw2();
     //hPERPBkg->Rebin(rebin);
     hPERPBkg->Scale(normC);
     TH1D* hPERPSub = new TH1D("hPERPSub", "Phi_Proton_pi_cutBE_sub", 100, minx, maxx);
     hPERPSub->Sumw2();
     hPERPSub->Add(hPERP, hPERPBkg, 1, -1); 
     hPERPSub->Draw();


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

     hPERPSub->Rebin(rebin);
     hPARASub->Rebin(rebin);

     Asym_m2g_bin[m-1][i][j][jj] = drawBeamAsymmetry_bin(hPERPSub, hPARASub, bin, norm, minx, maxx, color[2], figName, histName_h);

     Err_m2g_bin[m-1][i][j][jj] = Errtmp;
     
     //cout<<"x[j]="<<x[j]<<"; y[j]="<<y[j]<<"; ey[j]="<<ey[j]<<"; nt="<<nt<<endl;

     delete hPERPSub;
     delete hPARASub;
     fPERP->Close();
     fPARA->Close();

     }
    }
   }
  }

/*
  cout<<"n="<<nt-1<<endl;
  cout<<"x[n]={";
  for(int i=1; i<nt-2; i++)
    {
     cout<<x[0][0][0][0][i]<<", ";
    }
  cout<<x[0][0][0][0][nt-2]<<"}"<<endl;

  cout<<"y[n]={";
  for(int i=1; i<nt-2; i++)
    {
     cout<<Asym_tbin[0][0][0][0][i]<<", ";
    }
  cout<<Asym_tbin[0][0][0][0][nt-2]<<"}"<<endl;

  cout<<"ey[n]={";
  for(int i=1; i<nt-2; i++)
    {
     cout<<Err_tbin[0][0][0][0][i]<<", ";
    }
  cout<<Err_tbin[0][0][0][0][nt-2]<<"}"<<endl;
*/

  drawbinnedPSigma3_sys(x, Asym_m2g_bin, Err_m2g_bin, nt-1, nrebin, nsysSc, ncutsRg, nPS, "figPhiP/m2gbin/PSigma");
 
  return;
}

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 + ".pdf");
  c->Print(figName + ".ps");
  
  delete c;
  delete leg;  

  return;
}

double drawomega(TString histName, TString figName)
{
  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);
  
    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[m-1][i][j][1][2];
  ////double PSbkg = Asymomega; 
  ////double ff = fraction;
  //double PSsig = (PStot-ff*PSbkg)/(1-ff);
  //cout<<"PSsig(eta not corrected)="<<PStot<<endl;
  //cout<<"PSsig(eta corrected)="<<PSsig<<endl;

  fdata->Close();
  finmc->Close();
  fomega->Close();

  return fraction;
}


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->FixParameter(1, -0.07774);
  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->FixParameter(2, -0.07774);
  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 drawAsymSysErrOme() {

  cout<<"beginning"<<endl;

  //int rebin = 5;
  //int bin   = 100/rebin;
  Double_t normC = 1.0/6.0;
  double norm = 1.05;  //for ana & systematic err 1.1656
  //double norm = 1.12581;  //for ana & systematic err 1.1656
  //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 hist2DName_head = "Phi_Proton_vs_IM2g_cutBE_sys_"; // add any TH1D histogram name
  //TString hist2DName_head = "Phi_Proton_vs_t_sys_"; // add any TH1D histogram name
  TString tailfileName = ".root";

  TString figsName_head = "./figPhiP/standard/asym_";

  //const int nsysSc=9;  //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_", "dEdx_"}; 

  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};	

  TString figsDir = "./figsomega/"; 
  TString histName_ome_head = "S_IM_2gamma_sys_";
  TString figName_ome_head  = figsDir + "Im_2g_omega";

  cout<<__LINE__<<endl;
  drawomegaAsymSysErr(norm);

  cout<<__LINE__<<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++)
	 {
         Asymbkg[m-1][i][j][k][l] = 0.0;
         frac[m-1][i][j][k][l] = 0.0;
         if(k==1&&l==2)
            {
            Asymbkg[m-1][i][j][k][l] = Asym_m2g_bin[m-1][i][j][5];
            }
	 }
      }   
     }   
    }
   }
  

  //for(int rebin = 1; rebin<11; rebin++)
  for(int m = 1; m<nrebin+1; m++)
   {
   //if(m>2)  continue;
   int rebin = rebin_i[m-1];
   int bin   = 100/rebin;
   for(int i=0; i<nsysSc; i++)
    {
    //if(i>2)  continue;
   //for(int i=0; i<1; 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;
      TString histName_ome = histName_ome_head + figsName_mid;
      TString figName_ome = figName_ome_head + "rebin_" + title_rebin[m-1] + titleSc[i] + titleRg[j] + titlePS[k] + figsName_sec[2];

      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;

      if(k==1&&l==2)
         {
         cout<<"histName_ome= "<<histName_ome<<";   figName_ome = "<<figName_ome<<endl;
         cout<<drawomega(histName_ome, figName_ome)<<endl;
         frac[m-1][i][j][k][l] = drawomega(histName_ome, figName_ome);
         }

         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);

         Asymsig[m-1][i][j][k][l] = Asym[m-1][i][j][k][l];

         if((l == 2)&&(k == 1))
           {
           cout<<"Before correction:  "<<endl;
           cout<<"Asym for eta before omega correction for "<<i<<", "<<j<<", "<<k<<" = "<<Asym[m-1][i][j][k][l]<<endl;
           cout<<"Asym for eta after omega correction for "<<i<<", "<<j<<", "<<k<<" = "<<Asymsig[m-1][i][j][k][l]<<endl;
           }
 

         Err[m-1][i][j][k][l] = Errtmp;       

         delete hPERPSub;
         delete hPARASub;

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

         }

        double PSbkg =0.0;

        if(k==1)
           {
           PSbkg = Asymbkg[m-1][i][j][k][2];
           //PSbkg = -0.1736;
           Asymsig[m-1][i][j][k][2] = (Asym[m-1][i][j][k][2] - frac[m-1][i][j][k][2]*PSbkg)/(1-frac[m-1][i][j][k][2]); 
           //double PSsig = (PStot-ff*PSbkg)/(1-ff);
           cout<<"After correction:  "<<endl;
           cout<<"Asym for eta before omega correction for "<<i<<", "<<j<<", "<<k<<" = "<<Asym[m-1][i][j][k][2]<<endl;
           cout<<"Asym for eta after omega correction for "<<i<<", "<<j<<", "<<k<<" = "<<Asymsig[m-1][i][j][k][2]<<endl;
           }
        }
      }
     } 
    }

/*
  //for(int i=0; i<1; i++)
  // {
  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;
         }
        }
       }
     } 
    }
*/

  double min_Asym_pi0_50bins[nsysSc], max_Asym_pi0_50bins[nsysSc], uncertainty_pi0_50bins[nsysSc],  min_Asym_eta_50bins[nsysSc], max_Asym_eta_50bins[nsysSc], uncertainty_eta_50bins[nsysSc];
  double min_Asym_pi0_50bins_all, max_Asym_pi0_50bins_all, uncertainty_pi0_50bins_all, min_Asym_eta_50bins_all, max_Asym_eta_50bins_all, uncertainty_eta_50bins_all;

  for(int i=0; i<nsysSc; i++)
   {
    min_Asym_pi0_50bins[i] = Asym[1][0][0][0][2];
    max_Asym_pi0_50bins[i] = Asym[1][0][0][0][2];
    min_Asym_eta_50bins[i] = Asymsig[1][0][0][1][2];
    max_Asym_eta_50bins[i] = Asymsig[1][0][0][1][2];
   }

  cout<<"Beam Asymmetry for pi0 (rebin =2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {  
    if(Asym[1][i][j][0][2]<min_Asym_pi0_50bins[i])
        {
	min_Asym_pi0_50bins[i] = Asym[1][i][j][0][2];
        }
    if(Asym[1][i][j][0][2]>max_Asym_pi0_50bins[i])
        {
	max_Asym_pi0_50bins[i] = Asym[1][i][j][0][2];
        }

    cout<<Asym[1][i][j][0][2]<<" +/- "<<Err[1][i][j][0][2]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

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

  cout<<"Beam Asymmetry for eta after correction from omega background (rebin =2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {   
    if(Asymsig[1][i][j][1][2]<min_Asym_eta_50bins[i])
        {
	min_Asym_eta_50bins[i] = Asymsig[1][i][j][1][2];
        }
    if(Asymsig[1][i][j][1][2]>max_Asym_eta_50bins[i])
        {
	max_Asym_eta_50bins[i] = Asymsig[1][i][j][1][2];
        }

    cout<<Asymsig[1][i][j][1][2]<<" +/- "<<Err[1][i][j][1][2]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
    
     }
   }

  min_Asym_pi0_50bins_all = min_Asym_pi0_50bins[0];
  min_Asym_eta_50bins_all = min_Asym_eta_50bins[0];
  max_Asym_pi0_50bins_all = max_Asym_pi0_50bins[0];
  max_Asym_eta_50bins_all = max_Asym_eta_50bins[0];
  for(int i=0; i<nsysSc; i++)
     {
     uncertainty_pi0_50bins[i] = (max_Asym_pi0_50bins[i] - min_Asym_pi0_50bins[i])/Asym[1][0][0][0][2];
     uncertainty_eta_50bins[i] = (max_Asym_eta_50bins[i] - min_Asym_eta_50bins[i])/Asym[1][0][0][1][2];
     cout<<"min_Asym_pi0_50bins for "<<titleSc[i]<<" = "<<min_Asym_pi0_50bins[i]<<endl;
     cout<<"max_Asym_pi0_50bins for "<<titleSc[i]<<" = "<<max_Asym_pi0_50bins[i]<<endl;
     cout<<"min_Asym_eta_50bins for "<<titleSc[i]<<" = "<<min_Asym_eta_50bins[i]<<endl;
     cout<<"max_Asym_eta_50bins for "<<titleSc[i]<<" = "<<max_Asym_eta_50bins[i]<<endl;
     cout<<"uncertainty_pi0_50bins for "<<titleSc[i]<<" = "<<uncertainty_pi0_50bins[i]<<endl;
     cout<<"uncertainty_eta_50bins for "<<titleSc[i]<<" = "<<uncertainty_eta_50bins[i]<<endl;
     if(min_Asym_pi0_50bins[i]<min_Asym_pi0_50bins_all)
	{
	min_Asym_pi0_50bins_all = min_Asym_pi0_50bins[i];
	}
     if(min_Asym_eta_50bins[i]<min_Asym_eta_50bins_all)
	{
	min_Asym_eta_50bins_all = min_Asym_eta_50bins[i];
	}
     if(max_Asym_pi0_50bins[i]>max_Asym_pi0_50bins_all)
	{
	max_Asym_pi0_50bins_all = max_Asym_pi0_50bins[i];
	}
     if(max_Asym_eta_50bins[i]>max_Asym_eta_50bins_all)
	{
	max_Asym_eta_50bins_all = max_Asym_eta_50bins[i];
	}
     }
  uncertainty_pi0_50bins_all = (max_Asym_pi0_50bins_all - min_Asym_pi0_50bins_all)/Asym[1][0][0][0][2];
  uncertainty_eta_50bins_all = (max_Asym_eta_50bins_all - min_Asym_eta_50bins_all)/Asym[1][0][0][1][2];
  cout<<"uncertainty_pi0_50bins for all"<<" = "<<uncertainty_pi0_50bins_all<<endl;
  cout<<"uncertainty_eta_50bins for all"<<" = "<<uncertainty_eta_50bins_all<<endl;

  double Pol_Beam[3] = {0.3145, 0.4180, 0.4005};
  double Err_Pol[3]  = {0.0212, 0.0100, 0.0103};

  double Pol_Para[3] = {0.3129, 0.4321, 0.4103};
  double Err_Para[3] = {0.0373, 0.0153, 0.0170};

  double Pol_Perp[3] = {0.3161, 0.4039, 0.3907}; 
  double Err_Perp[3] = {0.0199, 0.0128, 0.0118};

  for(int m = 1; m<nrebin+1; m++)
   {
   for(int i=0; i<nsysSc; i++)
    {
    for(int j=0; j<ncutsRg; j++)
     {   
     Sigma_pi0[m-1][i][j] = Asym[m-1][i][j][0][2]/Pol_Beam[2];
     Sigma_eta[m-1][i][j] = Asymsig[m-1][i][j][1][2]/Pol_Beam[2];
     Err_Sigma_pi0[m-1][i][j] = Err[m-1][i][j][0][2]/Pol_Beam[2];
     Err_Sigma_eta[m-1][i][j] = Err[m-1][i][j][1][2]/Pol_Beam[2];
     if(i==2)
        {
        Sigma_pi0[m-1][i][j] = Asym[m-1][i][j][0][2]/Pol_Beam[j];
        Sigma_eta[m-1][i][j] = Asymsig[m-1][i][j][1][2]/Pol_Beam[j];
        Err_Sigma_pi0[m-1][i][j] = Err[m-1][i][j][0][2]/Pol_Beam[j];
        Err_Sigma_eta[m-1][i][j] = Err[m-1][i][j][1][2]/Pol_Beam[j];
        }
     }
    }
   }

       

  cout<<"Sigma Beam Asymmetry for pi0 (rebin = 2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {  
    cout<<Sigma_pi0[1][i][j]<<" +/- "<<Err_Sigma_pi0[1][i][j]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

  cout<<"Sigma Beam Asymmetry for eta (rebin = 2) corrected from omega"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {  
    cout<<Sigma_eta[1][i][j]<<" +/- "<<Err_Sigma_eta[1][i][j]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

  double min_Sigma_pi0_50bins[nsysSc], max_Sigma_pi0_50bins[nsysSc], uncert_Sigma_pi0_50bins[nsysSc],  min_Sigma_eta_50bins[nsysSc], max_Sigma_eta_50bins[nsysSc], uncert_Sigma_eta_50bins[nsysSc];
  double min_Sigma_pi0_50bins_all, max_Sigma_pi0_50bins_all, uncert_Sigma_pi0_50bins_all, min_Sigma_eta_50bins_all, max_Sigma_eta_50bins_all, uncert_Sigma_eta_50bins_all;

  for(int i=0; i<nsysSc; i++)
   {
    min_Sigma_pi0_50bins[i] = Sigma_pi0[1][0][0];
    max_Sigma_pi0_50bins[i] = Sigma_pi0[1][0][0];
    min_Sigma_eta_50bins[i] = Sigma_eta[1][0][0];
    max_Sigma_eta_50bins[i] = Sigma_eta[1][0][0];
   }

  //calculate the min and max for "Beam Asymmetry Sigma for pi0/eta (rebin =2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {  
    if(Sigma_pi0[1][i][j]<min_Sigma_pi0_50bins[i])
        {
	min_Sigma_pi0_50bins[i] = Sigma_pi0[1][i][j];
        }
    if(Sigma_pi0[1][i][j]>max_Sigma_pi0_50bins[i])
        {
	max_Sigma_pi0_50bins[i] = Sigma_pi0[1][i][j];
        }
    if(Sigma_eta[1][i][j]<min_Sigma_eta_50bins[i])
        {
	min_Sigma_eta_50bins[i] = Sigma_eta[1][i][j];
        }
    if(Sigma_eta[1][i][j]>max_Sigma_eta_50bins[i])
        {
	max_Sigma_eta_50bins[i] = Sigma_eta[1][i][j];
        }
     }
   }

  min_Sigma_pi0_50bins_all = min_Sigma_pi0_50bins[0];
  min_Sigma_eta_50bins_all = min_Sigma_eta_50bins[0];
  max_Sigma_pi0_50bins_all = max_Sigma_pi0_50bins[0];
  max_Sigma_eta_50bins_all = max_Sigma_eta_50bins[0];
  for(int i=0; i<nsysSc; i++)
     {
     uncert_Sigma_pi0_50bins[i] = (max_Sigma_pi0_50bins[i] - min_Sigma_pi0_50bins[i])/Sigma_pi0[1][0][0];
     uncert_Sigma_eta_50bins[i] = (max_Sigma_eta_50bins[i] - min_Sigma_eta_50bins[i])/Sigma_eta[1][0][0];
     cout<<"min_Sigma_pi0_50bins for "<<titleSc[i]<<" = "<<min_Sigma_pi0_50bins[i]<<endl;
     cout<<"max_Sigma_pi0_50bins for "<<titleSc[i]<<" = "<<max_Sigma_pi0_50bins[i]<<endl;
     cout<<"min_Sigma_eta_50bins for "<<titleSc[i]<<" = "<<min_Sigma_eta_50bins[i]<<endl;
     cout<<"max_Sigma_eta_50bins for "<<titleSc[i]<<" = "<<max_Sigma_eta_50bins[i]<<endl;
     cout<<"uncert_Sigma_pi0_50bins for "<<titleSc[i]<<" = "<<uncert_Sigma_pi0_50bins[i]<<endl;
     cout<<"uncert_Sigma_eta_50bins for "<<titleSc[i]<<" = "<<uncert_Sigma_eta_50bins[i]<<endl;
     if(min_Sigma_pi0_50bins[i]<min_Sigma_pi0_50bins_all)
	{
	min_Sigma_pi0_50bins_all = min_Sigma_pi0_50bins[i];
	}
     if(min_Sigma_eta_50bins[i]<min_Sigma_eta_50bins_all)
	{
	min_Sigma_eta_50bins_all = min_Sigma_eta_50bins[i];
	}
     if(max_Sigma_pi0_50bins[i]>max_Sigma_pi0_50bins_all)
	{
	max_Sigma_pi0_50bins_all = max_Sigma_pi0_50bins[i];
	}
     if(max_Sigma_eta_50bins[i]>max_Sigma_eta_50bins_all)
	{
	max_Sigma_eta_50bins_all = max_Sigma_eta_50bins[i];
	}
     }
  uncert_Sigma_pi0_50bins_all = (max_Sigma_pi0_50bins_all - min_Sigma_pi0_50bins_all)/Sigma_pi0[1][0][0];
  uncert_Sigma_eta_50bins_all = (max_Sigma_eta_50bins_all - min_Sigma_eta_50bins_all)/Sigma_eta[1][0][0];
  cout<<"uncert_Sigma_pi0_50bins for all"<<" = "<<uncert_Sigma_pi0_50bins_all<<endl;
  cout<<"uncert_Sigma_eta_50bins for all"<<" = "<<uncert_Sigma_eta_50bins_all<<endl;
  

  cout<<"Beam Asymmetry for omega (rebin =2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {   
    cout<<Asym_m2g_bin[1][i][j][5]<<" +/- "<<Err_m2g_bin[1][i][j][5]<<";      ";
    if(j==2)   cout<<titleSc[i]<<endl;
     }
   }

  cout<<"fraction for omega (rebin =2)"<<endl;
  for(int i=0; i<nsysSc; i++)
   {
   for(int j=0; j<ncutsRg; j++)
    {   
    cout<<frac[1][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 before correction from omega background"<<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<<"Beam Asymmetry for eta with rebin=1, 2, 4, 5, 10 for standard cut condition after correction from omega background"<<endl;
   for(int m = 1; m<nrebin+1; m++)
     {
     cout<<Asymsig[m-1][0][0][1][2]<<" +/- "<<Err[m-1][0][0][1][2]<<";      "<<endl;
     }

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

   cout<<"fraction for omega with rebin=1, 2, 4, 5, 10 for standard cut condition"<<endl;
   for(int m = 1; m<nrebin+1; m++)
     {
     cout<<frac[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++)
      {   
         cout<<hN[m-1][i][j][0][2]<<endl;
         cout<<"NPERP_a="<<NPERP_a[m-1][i][j][0][2]<<";      ";
         cout<<"NPARA_a="<<NPARA_a[m-1][i][j][0][2]<<";      ";
         cout<<"norm_a="<<Norm_a[m-1][i][j][0][2]<<";"<<endl;
      }
     } 
    }


  return;
}