#include <vector>

void peak_study3_plot(Double_t pedestal=100)
{


   // Plot pedestal-subtraced peaks produced by peak_study3.C

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15); 

gStyle->SetLabelSize(0.05,"XYZ");
gStyle->SetTitleSize(0.05,"XYZ");
//gStyle->SetTitleOffset(1.5,"Y");
gStyle->SetNdivisions(505,"XYZ");


// gStyle->SetFillColor(0);
//
   char string[256];


  const int npts=6;

  Double_t npeds[npts]={1,2,4,8,16,32};

  Double_t peak_int[npts]={-0.05, 0.28, 0.40, 0.46, 0.48, 0.50};
  Double_t sigma_int[npts]={1.25, 0.95, 0.72, 0.57, 0.44, 0.33};
  Double_t peak_float[npts]={-0.05, 0.03, 0.02, 0.01, 0.00};
  Double_t sigma_float[npts]={1.25, 0.92, 0.65, 0.52, 0.42, 0.35};

  /*Double_t peak_true[npts]={-0.49, -0.49, -0.49, -0.51, -0.50};
  Double_t sigma_true[npts]={0.30, 0.27, 0.26,  0.30, 0.28};*/

   TGraph *gpeak_int = new TGraph(npts,npeds,peak_int);
   TGraph *gsigma_int = new TGraph(npts,npeds,sigma_int);
   TGraph *gpeak_float = new TGraph(npts,npeds,peak_float);
   TGraph *gsigma_float = new TGraph(npts,npeds,sigma_float);

   TCanvas *c1 = new TCanvas("c1", "c1",200,10,1000,1000);
   c1->Divide(1,2);
   c1->SetGridx();
   c1->SetGridy();
   // c1->SetLogx();

   // ouput legend

   TLegend *leg = new TLegend(0.55,0.7,0.85,0.9);
   Double_t xmin=0;
   Double_t xmax=35;

   c1->cd(2);
   gsigma_int->SetTitle("Monte Carlo: Pedestal-subtracted Peak, Ped #sigma=1.2 counts");
   gsigma_int->GetYaxis()->SetTitle("#sigma (counts)");
   gsigma_int->GetXaxis()->SetTitle("Number of samples in pedestal average");
   gsigma_int->GetYaxis()->SetRangeUser(0,1.5);
   gsigma_int->GetXaxis()->SetRangeUser(xmin,xmax);
   gsigma_int->SetMarkerStyle(20);
   gsigma_int->SetMarkerColor(2);
   // gsigma_int->SetMarkerSize(0.5);
   gsigma_int->Draw("Ap");

   gsigma_float->SetMarkerStyle(20);
   gsigma_float->SetMarkerColor(4);
   // gsigma_float->SetMarkerSize(0.5);
   gsigma_float->Draw("samep");

  TLine *line = new TLine(xmin,1/sqrt(12),xmax,1/sqrt(12));
   line->Draw();

   sprintf (string,"#sigma=1/#sqrt{12}\n");
   printf("string=%s",string);
   t1 = new TLatex(2,0.32,string);
   t1->Draw();


   leg->AddEntry(gsigma_int,"Integer pedestal","p");
   leg->AddEntry(gsigma_float,"Floating pedestal","p");
   leg->Draw();

   c1->cd(1);
   gpeak_int->SetTitle("Monte Carlo: Pedestal-subtracted Peak, Ped #sigma=1.2 counts");
   gpeak_int->GetYaxis()->SetTitle("Computed Peak - Generated (counts)");
   gpeak_int->GetXaxis()->SetTitle("Number of samples in pedestal average");
   gpeak_int->GetYaxis()->SetRangeUser(-1,1.5);
   gpeak_int->GetXaxis()->SetRangeUser(xmin,xmax);
   gpeak_int->SetMarkerStyle(20);
   gpeak_int->SetMarkerColor(2);
   // gpeak_int->SetMarkerSize(0.5);
   gpeak_int->Draw("Ap");

  TLine *line = new TLine(xmin,0,xmax,0);
   line->Draw();

   gpeak_float->SetMarkerStyle(20);
   gpeak_float->SetMarkerColor(4);
   // gpeak_float->SetMarkerSize(0.5);
   gpeak_float->Draw("samep");

   TLegend *leg = new TLegend(0.55,0.7,0.85,0.9);

   leg->AddEntry(gpeak_int,"Integer pedestal","p");
   leg->AddEntry(gpeak_float,"Floating pedestal","p");
   leg->Draw();



 


   sprintf (string,"peak_study3_plot.pdf");
   c1->SaveAs(string);



}
Double_t pulse_func (Double_t *x, Double_t *par) {

   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t gain=par[2];
   Double_t t0=par[3];
   Double_t pedestal=par[4];
   Double_t x1=x[0]-t0;

   char string[256];

   Double_t func;
   Double_t amplitude;

   if (x1 < -20) {
           func = pedestal;
           }
   else if {

         Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
         Double_t arg2 = -mu*x1 +(mu*sigma)*(mu*sigma)/2;
         amplitude = (mu/2)* exp(arg2) *  TMath::Erfc(arg);
          func = gain * amplitude + pedestal;
  }
   
   /*sprintf (string,"x1=%f amplitude=%f mu=%f sigma=%f t0=%f pedestal=%f, func=%f\n",x1,amplitude,mu,sigma,pedestal,func);
   printf ("string=%s",string);*/

   return func;
}