#include "/home/zhaobo/top_dir/include/histerror.C"

#define deg2rad 180/3.1416
#define N 4
#define lowx1 0.1
#define highx1 0.6

double background(double *x, double *par)
{
 double pol = 0;
 for(int i=0; i<N; i++)
  pol = pol + par[i]*pow(x[0], i);

 return pol;
}


double gausfun(double *x, double *par)
{
 if(par[2] != 0) double arg1 = (x[0] - par[1])/par[2];
 return par[0]*TMath::Exp(-0.5*arg1*arg1);
}

double fitf(double *x, double *par)
{
 return background(x, par) + gausfun(x, &par[N]);
}


double fitf1(double *x, double *par)
{
 return (par[0]*sin(x[0]/deg2rad))/(1 + par[1]*cos(x[0]/deg2rad) + par[2]*cos(2*x[0]/deg2rad));
}

double fitf2(double *x, double *par)
{
 return par[0]*sin(x[0]/deg2rad) + par[1]*sin(2*x[0]/deg2rad);
}


void etaplotbsa_MM2()
{
 gStyle->SetOptFit(1);
 gStyle->SetPalette(1);
 gStyle->SetStatBorderSize(2);
 gStyle->SetStatW(0.14);
 gStyle->SetStatH(0.19);

 TLatex title;
 title.SetTextColor(4);
 title.SetTextFont(52);
 title.SetTextSize(0.04);
 title.SetNDC();


 TFile *infile = new TFile("allhists_MM2.root", "READ");

 TH1F *heta_MM2_phi_sign[12][2]; 
 TH1F *heta_wtMM2_phi_sign[12][2];

  for(int l=0; l<12; l++)
   for(int n=0; n<2; n++)
   {
    heta_MM2_phi_sign[l][n] = (TH1F*)infile->Get(Form("heta_MM2_phi%d_sign%d", l, n));
    heta_wtMM2_phi_sign[l][n] = (TH1F*)infile->Get(Form("heta_wtMM2_phi%d_sign%d", l, n));
   }

 TFile *outf = new TFile("etafinalbsa_MM2.root", "RECREATE");
 TH1F *hbsa_eta;
 TH1F *hist_eta[2], *histwt_eta[2];

  for(int n=0; n<2; n++)
  {
   hist_eta[n] = new TH1F(Form("hist_eta_sign%d", n), "each #phi bin", 12, 0, 360);
   histwt_eta[n] = new TH1F(Form("histwt_eta_sign%d", n), "each #phi bin", 12, 0, 360);
  }

// TF1 *fgaus = new TF1("fgaus", gausfun, 0.1, 0.6, 3);

// fgaus->SetParameters(50, 0.3, 0.05);
// fgaus->SetParLimits(2, 0, 1); 
 double par[N+3];

 TF1 *backFcn1 = new TF1("backFcn1", background, lowx1, highx1, N);
 backFcn1->SetLineColor(kBlue);
 backFcn1->SetLineWidth(1);

 TF1 *gausFcn1 = new TF1("gausFcn1", gausfun, lowx1, highx1, 3);
 gausFcn1->SetLineColor(kGreen);
 gausFcn1->SetLineWidth(1);

 TF1 *func1 = new TF1("func1", fitf, lowx1, highx1, N+3);
 func1->SetLineColor(kRed);
 func1->SetLineWidth(1);


 func1->SetParameters(0, 0, 0, 0, 3e1, 0.3, 0.1);
 func1->SetParLimits(4, 0, 200);
 func1->SetParLimits(5, 0.2, 0.4);
 func1->SetParLimits(6, 0, 1);


 TCanvas *c1 = new TCanvas("c1", "c1", 800, 600);
 c1->Divide(3, 4);

 c1->Print("etafitMM2.ps[");
  for(int n=0; n<1; n++)
  {
   for(int l=0; l<12; l++)
   {
    c1->cd(l+1);
    heta_MM2_phi_sign[l][n]->Fit("func1", "EMR");
    func1->GetParameters(par);
    backFcn1->SetParameters(par);
    backFcn1->Draw("same");
    gausFcn1->SetParameters(&par[N]);
    gausFcn1->Draw("same");
    gPad->Update();
//    c1->Update();
//    c1->Modified();
    hist_eta[n]->SetBinContent(l+1, par[N]);
    hist_eta[n]->SetBinError(l+1, func1->GetParError(N));
   }   
   c1->Print("etafitMM2.ps");
  } 
  c1->Print("etafitMM2.ps]");

/*
  c1->Print("etafitMM2wt.ps[");
  for(int n=0; n<2; n++)
  {
   for(int l=0; l<12; l++)
   {
    c1->cd(l+1);
    heta_wtMM2_phi_sign[l][n]->Fit("func1", "EMR");
    func1->GetParameters(par);
    histwt_eta[n]->SetBinContent(l+1, par[N]);
    histwt_eta[n]->SetBinError(l+1, func1->GetParError(N));
   }
   c1->Print("etafitMM2wt.ps");
  }
  c1->Print("etafitMM2wt.ps]");
*/

 outf->Write();

}