#include <vector>

void peak_study3(Double_t pedestal=100)
{


   // Determine optimal peak - pedestal algorithm comparing floating and integer operations

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;
   Int_t nevents=1000;

    Double_t ydelta=10;

    Double_t peak0=249;
    Double_t Dpeak=0.5;
    // Double_t Dpeak=1.5;
    Double_t pedestal=100;
    Int_t nped=16;

   TH1D* hp_true = new TH1D("hp_true","Generated peak",nbins, 250-ydelta,250+ydelta);
   TH1D* hped_true = new TH1D("hped_true","Generated pedestal",nbins, 100-ydelta,100+ydelta);
   TH1D* hpulse = new TH1D("hpulse","Generated pulse",nbins, 150-ydelta,150+ydelta);

   TH1D* hp_data = new TH1D("hp_data","Data peak",nbins, 250-ydelta,250+ydelta);
   TH1D* hped_firm = new TH1D("hped_firm","Firmware pedestal",nbins, 100-ydelta,100+ydelta);
   TH1D* hpulse_firm = new TH1D("hpulse_firm","Firmware pulse",nbins, 150-ydelta,150+ydelta);

   TH1D* hpulse_float = new TH1D("hpulse_float","Pulse-float pedestal",nbins, 150-ydelta,150+ydelta);

    Double_t ydelta=5;
   TH1D* hdiff_firm = new TH1D("hdiff_firm","Pulse firm - Pulse true",nbins, -ydelta,ydelta);
   TH1D* hdiff_float = new TH1D("hdiff_float","Pulse float - Pulse true",nbins, -ydelta,ydelta);

   Double_t ped_sigma=1.2;
   TRandom3 *r1 = new TRandom3();
   UInt_t  iseed = 23424111513;
   r1->SetSeed(iseed);

   Int_t Nbits=5;
   Int_t npeds = pow(2,Nbits);
   cout << " npeds=" << npeds << endl;


   // loop over nevents 

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

     // randomize peak values
      Double_t p_true = peak0 - Dpeak + 2*Dpeak*r1->Uniform() + r1->Gaus(0,ped_sigma);
      Int_t p_data = p_true;

      // here we need to add loop to obtain pedestal sum

     Double_t ped_sum=0;
     Int_t iped_sum=0;

       for (Int_t jj=0; jj<npeds;jj++) {
           Double_t ped = r1->Gaus(pedestal,ped_sigma);
           ped_sum += ped;
           iped_sum += ped;
           }
           
      Double_t ped_true = ped_sum/npeds;
      Int_t ped_firm = iped_sum >> Nbits;

      // cout << " j=" << j << " ped_true=" << ped_true << " ped_firm=" << ped_firm << endl;

      Double_t pulse = p_true - pedestal;
      Double_t ped_float = iped_sum/(float) npeds;

      Int_t pulse_firm = p_data - ped_firm;
      Double_t pulse_float = p_data - ped_float;

     // cout << " j=" << j << " p_true=" << p_true << " Dpeak=" << Dpeak << " ped_true=" << ped_true << " pulse=" << pulse << " p_data=" << p_data << " ped_firm=" << ped_firm << " ped_float=" << ped_float << " pulse_firm=" << pulse_firm << " pulse_float=" << pulse_float << endl;

      // fill histograms

     hp_true->Fill (p_true);
     hped_true->Fill(ped_true);
     hpulse->Fill(pulse);
   
     hp_data->Fill(p_data);
     hped_firm->Fill(ped_firm);
     hpulse_firm->Fill(pulse_firm);
  
     hpulse_float->Fill(pulse_float);

     hdiff_firm->Fill(pulse_firm - pulse);
     hdiff_float->Fill(pulse_float - pulse);
   }

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

   // ouput legend

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

   c1->cd(1);
   hp_true->SetXTitle("ADC counts");
   hp_true->SetYTitle("Events");
   hp_true->Draw();

   sprintf (string,"Gen peak=%.1f-%.1f\n",peak0-Dpeak,peak0+Dpeak);
   printf("string=%s",string);
   t1 = new TLatex(0.55,0.6,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Ped mean=%d\n",pedestal);
   printf("string=%s",string);
   t1 = new TLatex(0.55,0.56,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Ped width=%.1f\n",ped_sigma);
   printf("string=%s",string);
   t1 = new TLatex(0.55,0.52,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Nsamples=%d\n",npeds);
   printf("string=%s",string);
   t1 = new TLatex(0.55,0.48,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c1->cd(2);
   hped_true->SetXTitle("ADC counts");
   hped_true->SetYTitle("Events");
   hped_true->Draw();

   c1->cd(3);
   hpulse->SetXTitle("ADC counts");
   hpulse->SetYTitle("Events");
   hpulse->Draw();

   c1->cd(4);
   hp_data->SetXTitle("ADC counts");
   hp_data->SetYTitle("Events");
   hp_data->Draw();

   c1->cd(5);
   hped_firm->SetXTitle("ADC counts");
   hped_firm->SetYTitle("Events");
   hped_firm->Draw();

   c1->cd(6);
   hpulse_firm->SetXTitle("ADC counts");
   hpulse_firm->SetYTitle("Events");
   hpulse_firm->Draw();

   c1->cd(7);
   hpulse_float->SetXTitle("ADC counts");
   hpulse_float->SetYTitle("Events");
   hpulse_float->Draw();

   c1->cd(8);
   hdiff_firm->SetXTitle("ADC counts");
   hdiff_firm->SetYTitle("Events");
   hdiff_firm->Draw();

   sprintf (string,"Gen peak=%.1f-%.1f\n",peak0-Dpeak,peak0+Dpeak);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Ped mean=%d\n",pedestal);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.76,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Ped width=%.1f\n",ped_sigma);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.72,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Nsamples=%d\n",npeds);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.68,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c1->cd(9);
   hdiff_float->SetXTitle("ADC counts");
   hdiff_float->SetYTitle("Events");
   hdiff_float->Draw();

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

 


   sprintf (string,"peak_study3_%02d.pdf",npeds);
   c1->SaveAs(string);
   /*sprintf (string,"peak_study3_%02d.pdf)",npeds);
   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;
}