#include <vector>

void peak_study2(Double_t pedestal=100)
{

   // loop over generated pulses with variable heights and starting points relative to sampling and compute mean and sigmas for
  // various estimates of the peak

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


   // define histograms of peak estimates 

   Int_t nbins=100;

    Double_t ypmin=20;
    Double_t ypmax=30;

   TH1D* hptrue = new TH1D("hptrue","Generated peak",nbins, ypmin,ypmax);
   TH1D* hpsample = new TH1D("hpsample","Single sample max",nbins, ypmin,ypmax);
   TH1D* hp3sample = new TH1D("hp3sample","Average 3 max samples",nbins, ypmin,ypmax);
   TH1D* hpeak = new TH1D("hpeak","Parabolic peak estimate, all double",nbins, ypmin,ypmax);
   TH1D* hipeak = new TH1D("hipeak","Parabolic ipeak, double ratio",nbins, ypmin,ypmax);

   Double_t pedestal=2;
   Double_t ped_sigma=1.2;

   Double_t xpmin=0;
   Double_t xpmax=100;

   Double_t mu=0.0787;    // unis of 1/x
   Double_t sigma=0.865;     // units of x
   Double_t gain0=350;   
   Double_t t0=39.5;

   TF1 *pulse = new TF1("pulse_func",pulse_func,xpmin,xpmax,5);
   pulse->SetParNames("mu","sigma","gain","t0","pedestal");

   Double_t gain_min=gain0*1.00;
   Double_t gain_max=gain0*1.04;
   Int_t ngain = 10;
   Int_t nsample = 10;

   for (Int_t j=0;j<ngain;j++) {

   Double_t gain = gain_min + (gain_max-gain_min)*j/ngain;

   for (Int_t k=0; k<nsample;k++) {

   Double_t t = t0 + k/10.;   // taking fractions of one sample

    //cout << " gain=" << gain << "  t=" << t << endl;
   
   pulse->SetParameters(mu,sigma,gain,t,pedestal);

  // data from function

   const int npts=100;
  Double_t x[npts];
  Double_t y[npts];
  Int_t ix[npts];
  Int_t iy[npts];

 const int ntop=3;
 Double_t xtop[npts]={0,0,0};
 Double_t ytop[npts]={0,0,0};
 Int_t ixtop[npts]={0,0,0};
 Int_t iytop[npts]={0,0,0};

 Double_t tsample=1;   // 4 ns samples
 Double_t ypmax = 0;
 Double_t iypmax = 0;
 
  for (Int_t jj=1;jj<npts;jj++) {
         x[jj] = xpmin + jj*tsample;
          y[jj] = pulse->Eval(x[jj]);
          if ( ypmax < y[jj]) {
               xtop[0] = x[jj-1];
               xtop[1] = x[jj];
               xtop[2] = xpmin + (jj+1)*tsample;
               ytop[0] = y[jj-1];
               ytop[1] = y[jj];
               ytop[2] = pulse->Eval(xtop[2]);
               ypmax = y[jj];
              }

         // integerize

         ix[jj] = (int)x[jj];
         iy[jj] = (int)y[jj];
          if ( iypmax < iy[jj]) {
               ixtop[0] = ix[jj-1];
               ixtop[1] = ix[jj];
               ixtop[2] = xpmin + (jj+1)*tsample;
               iytop[0] = iy[jj-1];
               iytop[1] = iy[jj];
               iytop[2] = (int)pulse->Eval(ixtop[2]);
               iypmax = iy[jj];
              }
    }

  // Get maximum

  Double_t pmax = pulse->GetMaximum();

  // cout << " ytop[0] =" << ytop[0] << " ytop[1] =" << ytop[1] << " ytop[2] =" << ytop[2] << endl;
  // cout << " iytop[0] =" << iytop[0] << " iytop[1] =" << iytop[1] << " iytop[2] =" << iytop[2] << endl;

  // cout << "pmax=" << pmax << " ypmax=" << ypmax << " iypmax=" << iypmax << endl;
  

  TGraph *points = new TGraph(npts,ix,iy);

   // compute maximum


  Double_t x_peak = xtop[1] -0.5 * (ytop[2]-ytop[0])/(ytop[0]-2*ytop[1]+ytop[2]);
  Double_t y_peak = ytop[1] - (ytop[2]-ytop[0])*(ytop[2]-ytop[0])/(ytop[0]-2*ytop[1]+ytop[2])/8;

  Double_t ix_peak = ixtop[1] -0.5 * (iytop[2]-iytop[0])/(float)(iytop[0]-2*iytop[1]+iytop[2]);
  Double_t iy_peak = iytop[1] - (iytop[2]-iytop[0])*(iytop[2]-iytop[0])/(float)(iytop[0]-2*iytop[1]+iytop[2])/8;

  // cout << "x_peak=" << x_peak << " y_peak=" << y_peak << endl;
  // cout << "ix_peak=" << ix_peak << " iy_peak=" << iy_peak << endl;

  // Fill histograms

  hptrue->Fill(pmax);
  hpsample->Fill(iytop[1]);
  hp3sample->Fill( (iytop[0]+iytop[1]+iytop[2])/3.);
  hpeak->Fill(y_peak);
  hipeak->Fill(iy_peak);

   }  // end loop over samples


   }   // end loop over gains

  TMarker *peak = new TMarker(ix_peak, iy_peak,21);

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

   // ouput legend

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

   pulse->Draw();
   pulse->SetLineColor(1);
   points->SetMarkerStyle(20);
   points->SetMarkerColor(4);
   points->SetMarkerSize(0.7);
   points->Draw("psame");

   peak->SetMarkerStyle(21);
   peak->SetMarkerColor(2);
   peak->SetMarkerSize(1.2);
   peak->Draw("psame");

   leg->AddEntry(pulse,"Generated pulse","l");
   leg->AddEntry(points,"Samples","p");
   leg->AddEntry(peak,"Parabolic peak estimate","p");
   leg->Draw();

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

   c2->cd(1);
   hptrue->Draw();

   c2->cd(2);
   hpsample->Draw();

   c2->cd(3);
   hp3sample->Draw();
   
   c2->cd(4);
   hpeak->Draw();
   
   c2->cd(5);
   hipeak->Draw();

   c2->cd(6);
   sprintf (string,"For parabolic peak estimate:\n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.9,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"(x_{2},y_{2}) = max sample\n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"y_{p} = y_{2} - #frac{1}{8} (y_{3}-y_{1})^{2}/(y_{1}-2y_{2}+y_{3}) \n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.7,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"x_{p} = x_{2} - #frac{1}{2} (y_{3}-y_{1})/(y_{1}-2y_{2}+y_{3}) \n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.6,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"peak_study2.pdf(");
   c1->SaveAs(string);
   sprintf (string,"peak_study2.pdf)");
   c2->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;
}