#include <vector>

void pixel_size(void)
{

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->SetFillColor(0);
//
   char string[256];
   Int_t j,jj;

  Int_t nbins=100;

   // parameter definitions

   Double_t mu=0.0787;    // unis of 1/x
   Double_t sigma=0.865;     // units of x 
   Double_t t0=10;
   Double_t pedestal=100;
   Double_t ped_sigma=1.2; 
   Double_t Vover_nominal = 1.2; 
   Double_t Vover0 = 0.9;
   Double_t T0 = 18;
   Double_t T=10;
   Double_t gain =   0.75e6;                   // gain = 0.75x10^6 at 23 deg C and Vover=0.9 V.
   Double_t e = 1.6e-19;   // electric charge
   Double_t preamp_gain = 1.62*50;    // units of Ohms                // preamp gain 0.081mV/uA or x1.62 wrt to 50Ohm (Fernando 4/15/2015)
   Double_t fadc_scale = 2./4095.;      // 2 Volts for a 12 bit scale.  (V/count)
   Double_t time_sample = 4e-9;        // 4 ns time sample
   Double_t conversion = 0.046;         // Average MeV per integral count.
   Double_t norm_factor=1;               // relative normalization to play with gains

   Double_t pe_pulse = norm_factor*e * preamp_gain * gain / fadc_scale / time_sample; 
   Double_t norm = pe_pulse*Vover_nominal/Vover0;
   // Double_t norm=1;
 
   cout << " SiPM gain=" << gain << " , fadc gain=" << preamp_gain << " Ohms, fadc scale=" << fadc_scale << "V/count, time sample=" << time_sample << " ns, pe_pulse=" << pe_pulse << " counts" << " norm=" << norm << endl;

   Double_t xmin=0;
   Double_t xmax=100;

   Double_t xpmin=0;
   Double_t xpmax=100;

   Int_t npts=16;
          
   // define histogram to hold pulses

  TH1F *hpedestalw = new TH1F("hpedestalw","Pedestal window", nbins,xmin,xmax);
  TH1F *hpulsew = new TH1F("hpulsew","Pulse waveform", nbins,xmin,xmax);
  TH1F *hwindow= new TH1F("hwindow","Readout Window", nbins,xmin,xmax);
  TH1I *hwindowI = new TH1I("hwindowI","Readout Window Integer", nbins,xmin,xmax);
  
   // define pulse function: time units are FADC samples.

   TF1 *pulse = new TF1("pulse_func",pulse_func,xpmin,xpmax,5);
   pulse->SetParameters(mu,sigma,norm,t0,pedestal);
   pulse->SetParNames("mu","sigma","norm","t0","pedestal");
   TRandom1 *r = new TRandom1();      // randomize pedestals with FADC noise

        // define pulse function
 
         TF1 *pulse = new TF1(string,pulse_func,xpmin,xpmax,5);
         pulse->SetParameters(mu,sigma,norm,t0,pedestal);
         pulse->SetParNames("mu","sigma","norm","t0","pedestal");

   // create histogram windows with one pulse and electronic noise

  for (jj=0;jj<nbins;jj++) {
         Double_t t = xmin + (xmax-xmin)*(jj+0.5)/nbins;
         Double_t ran_ped = r->Gaus(0,ped_sigma);
          Double_t amp = pulse->Eval(t) - pedestal;
         Double_t fadc = amp + ran_ped;

             cout << " jj=" << jj << " t=" << t << " amp=" << amp << " ran_ped=" << ran_ped << " fadc=" << fadc << endl;
             hpedestalw->Fill(t,ran_ped);
             hpulsew->Fill(t,amp);
             hwindow->Fill(t,fadc);
             hwindowI->Fill(t,fadc+0.5);
         }
  
   TCanvas *c0 = new TCanvas("c0", "c0",200,10,700,700);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);

   c0->SetGridx();
   c0->SetGridy();

   c0->Divide(1,2);
   c0->cd(1);
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);

   hpulsew->SetTitle("");  
   hpulsew->GetXaxis()->SetLabelSize(0.05);
   hpulsew->GetXaxis()->SetTitleSize(0.05);
   hpulsew->GetYaxis()->SetLabelSize(0.05);
   hpulsew->GetYaxis()->SetTitleSize(0.05);
   hpulsew->GetYaxis()->SetTitleOffset(1.5);
   hpulsew->GetYaxis()->SetTitle("One pixel (FADC counts/sample)");
   hpulsew->GetXaxis()->SetTitle("FADC samples");
   hpulsew->GetXaxis()->SetNdivisions(505);
   hpulsew->SetLineColor(2);
   hpulsew->Draw();  

   sprintf (string,"Parameter Settings\n");
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.86,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"SiPM Gain=%g at Vover=%.1f V\n",gain,Vover0);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Preamp Gain=%.1f Ohms\n",preamp_gain);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"FADC scale=%.7f V/count\n",fadc_scale);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"Conversion=%.4f MeV/integral count\n",conversion);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.65,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"One pixel integral =%.2f counts\n",norm);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.61,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   sprintf (string,"One pixel integral =%.2f MeV\n",norm*conversion);
   printf("string=%s",string);
   t1 = new TLatex(0.4,0.57,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c0->cd(2);
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);

   hpedestalw->SetTitle("");
   hpedestalw->GetXaxis()->SetLabelSize(0.05);
   hpedestalw->GetXaxis()->SetTitleSize(0.05);
   hpedestalw->GetYaxis()->SetLabelSize(0.05);
   hpedestalw->GetYaxis()->SetTitleSize(0.05);
   hpedestalw->GetYaxis()->SetTitleOffset(1.5);
   hpedestalw->GetYaxis()->SetTitle("Random pedestal (counts/sample)");
   hpedestalw->GetXaxis()->SetTitle("FADC samples");
   hpedestalw->GetXaxis()->SetNdivisions(505);
   hpedestalw->SetLineColor(2);
   hpedestalw->Draw();      
  
   TCanvas *c1 = new TCanvas("c1", "c1",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();

   c1->Divide(1,2);
   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);

   hwindow->SetTitle("");
   hwindow->GetXaxis()->SetLabelSize(0.05);
   hwindow->GetXaxis()->SetTitleSize(0.05);
   hwindow->GetYaxis()->SetLabelSize(0.05);
   hwindow->GetYaxis()->SetTitleSize(0.05);
   hwindow->GetYaxis()->SetTitleOffset(1.5);
   hwindow->GetYaxis()->SetTitle("Float Counts");
   hwindow->GetXaxis()->SetTitle("FADC samples");
   hwindow->GetXaxis()->SetNdivisions(505);
   hwindow->SetLineColor(2);
   hwindow->Draw();   

   /*sprintf (string,"t0=%.1f cm, R=%.2f\n",xpos,Rup);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();*/

   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);

   hwindowI->SetTitle("");
   hwindowI->GetXaxis()->SetLabelSize(0.05);
   hwindowI->GetXaxis()->SetTitleSize(0.05);
   hwindowI->GetYaxis()->SetLabelSize(0.05);
   hwindowI->GetYaxis()->SetTitleSize(0.05);
   hwindowI->GetYaxis()->SetTitleOffset(1.5);
   hwindowI->GetYaxis()->SetTitle("Integer Counts");
   hwindowI->GetXaxis()->SetTitle("FADC samples");
   hwindowI->GetXaxis()->SetNdivisions(505);
   hwindowI->SetLineColor(2);
   hwindowI->Draw();       
  




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

   // save width vs rate histogram.

   sprintf (string,"pixel_size_c0.C");
   c0->SaveAs(string);
   sprintf (string,"pixel_size_c1.C");
   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,"pulse_func: 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;
}