#include <vector>

void pedestal_width(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;
   Double_t mu=0.0787;    // unis of 1/x
   Double_t sigma=0.865;     // units of x 
   Double_t gain=40;   
   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;
   Int_t layer=1;
   // Double_t gain =   (Vover_nominal/Vover0)*0.75e6;                   // gain = 0.75x10^6 at 23 deg C and Vover=0.9 V.
   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 norm_factor=1;                // relative to nominal to study pulse height dependence.
   Double_t pe_pulse = norm_factor*e * preamp_gain * gain / fadc_scale / time_sample; 
   Double_t norm = pe_pulse;
   // 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;

   Double_t fudge_dr=1.4;    // used to match measured pedestal width data
   Double_t dark_rate0=fudge_dr*layer*16*0.75e6*exp(0.075*(T-T0));       // 16 times rate per cell.  at Vover0; Eq 7 of NIMA698(2013)234.
   Double_t xpbegin=60;        // start sequence before actual window to include pulses before window start
   Double_t dark_mean0= dark_rate0*(xpmax-xpmin+xpbegin)*4e-9;     // average number of hits per window.

   // first fit data in NIMA 698 (2013) 234, fig 8. Read some values off the plot. Note that this section no longer used. 
  Int_t npts=5;
  Float_t vbr=69.5;
  Float_t dc_nim[5]={0,3.6,11.5,19.0,27.2};
  Float_t v_nim[5]={69.5,70.0,70.5,70.8,71.};
  for (j=0;j<npts;j++) {
          v_nim[j] -= vbr;      // compute voltage over breakdown
          if (j>0) dc_nim[j] /=v_nim[j];     // dark current needs to be divided by gain (proportional to overvoltage)  to get a rate. This should be linear vbr.
          cout << "j=" << j << " v_nim=" << v_nim[j] << " dc_nim=" << dc_nim[j] << endl;
          }

   TGraph *Delta_dr = new TGraph(npts,v_nim,dc_nim);
   Delta_dr->Fit("pol2");

   // get parameters
  
   TF1 *fit = Delta_dr->GetFunction("pol2");
   Float_t par0 = fit->GetParameter(0); 
   Float_t par1 = fit->GetParameter(1); 
   Float_t par2 = fit->GetParameter(2); 
   cout << " Fit parameters: par0=" << par0 << " par1=" << par1 << " par2=" << par2 << endl;

   Double_t Vover=1.2;
   // Double_t dark_rate = dark_rate0*(par0 + Vover*par1 + Vover*Vover*par2)/(par0 + Vover_nominal*par1 + Vover_nominal*Vover_nominal*par2);
   // Double_t dark_rate= dark_rate0*Vover/Vover_nominal;    // assume rate is proportional with voltage
   Double_t dark_rate= dark_rate0*Vover/Vover0;    // assume rate is proportional with voltage
   cout << "Dark rate = " << dark_rate << " at Vover=" << Vover << endl;

   TCanvas *cdr = new TCanvas("cdr", "cdr",200,10,700,700);
   cdr->SetBorderMode(0);
   cdr->SetFillColor(0);

   cdr->SetGridx();
   cdr->SetGridy();

   Delta_dr->SetTitle("");
   Delta_dr->GetXaxis()->SetLabelSize(0.05);
   Delta_dr->GetXaxis()->SetTitleSize(0.05);
   Delta_dr->GetYaxis()->SetLabelSize(0.05);
   Delta_dr->GetYaxis()->SetTitleSize(0.05);
   Delta_dr->GetYaxis()->SetTitleOffset(1.5);
   Delta_dr->GetYaxis()->SetTitle("Dark Rate (MHz)");
   Delta_dr->GetXaxis()->SetTitle("Over Voltage (V)");
   Delta_dr->GetXaxis()->SetNdivisions(505);
   Delta_dr->SetLineColor(2);
   Delta_dr->SetMarkerStyle(21);
   Delta_dr->Draw("AP");   

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

  TH1F *hhits= new TH1F("hhits","Number of dark hits/window", nbins,0,20);
  TH1F *hpulse= new TH1F("hpulse","Signal", nbins,0,5*ped_sigma);
  TH1F *hpedestal= new TH1F("hpedestal","Pedestal", nbins,-5*ped_sigma,5*ped_sigma);
  TH1F *hsamples= new TH1F("hsamples","All samples", nbins,pedestal-5*ped_sigma,pedestal+15*ped_sigma);

  
  Int_t kmax=11;
  TH2F *h2_ samples= new TH2F("h2_samples","Samples vs Vover", kmax,0,2.2,nbins,pedestal-5*ped_sigma,pedestal+15*ped_sigma);
  Double_t rmin=0;
  Double_t rmax=150;
  TH2F *h2_ rate= new TH2F("h2_rate","Samples vs rate", 20,rmin,rmax,nbins,pedestal-5*ped_sigma,pedestal+15*ped_sigma);

   // define function for sigma to rate

   TF1 *sqroot = new TF1("sqroot_func",sqroot_func,rmin,rmax,2);
   Double_t parA=1.2;
   Double_t parB=1.2e-2;
   sqroot->SetParameters(parA,parB);
   sqroot->SetParNames("parA","parB");

   // 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");

   vector <TF1*> pulses;

   TRandom1 *r = new TRandom1();      // randomize pedestals with FADC noise

   // pick actual number of hits in this window
   
   /*TCanvas *cdump = new TCanvas("cdump", "cdump",200,10,700,700);
   cdump->SetBorderMode(0);
   cdump->SetFillColor(0);
   cdump->Divide(2,3);*/
 
  for (Int_t k=0; k<kmax; k++) {       // loop over over voltage

  Double_t Vover= 0.2*k + 0.1;     // put value at center of bin
  norm = pe_pulse*Vover/Vover0;

  // Double_t dark_rate = dark_rate0*(par0 + Vover*par1 + Vover*Vover*par2)/(par0 + Vover_nominal*par1 + Vover_nominal*Vover_nominal*par2);
  Double_t dark_rate= dark_rate0*Vover/Vover0;    // assume rate is proportional with voltage
  Double_t dark_mean= dark_rate*(xpmax-xpmin+xpbegin)*4e-9;     // average number of hits per window.

  Int_t ncombos=100;

  for (Int_t kk=0;kk<ncombos; kk++) {       // take ncombos combinations per setting

  Int_t nhits = r->Poisson(dark_mean);
  // cout << "Dark_rate=" << dark_rate/1e6 << " MHz, window=" << 4*(xpmax-xpmin) << " ns, dark_mean=" << dark_mean << " nhits=" << nhits << " kk=" << kk << " kmax=" << kmax << endl;

   for(unsigned int i=0; i<pulses.size(); i++) delete pulses[i];
   pulses.clear();
   for (jj=0; jj<nhits;jj++) {
          sprintf (string,"clone%d",jj);
          Double_t tpulse = -xpbegin + r->Uniform()*(xpmax-xpmin+xpbegin);  
         TF1 *pulse = new TF1(string,pulse_func,xpmin,xpmax,5);
         pulse->SetParameters(mu,sigma,norm,tpulse,pedestal);
         // pulse->SetParNames("mu","sigma","norm","t0","pedestal");

          /*// cdump->cd(jj+1);
          // TF1 * pulse2 = pulse->DrawCopy();
          sprintf (string,"clone%d",jj);
          TF1 * pulse2 = (TF1*)pulse->Clone(string);
          pulse2->SetParameters(mu,sigma,norm,tpulse,pedestal);
          // cout << " jj=" << jj << " tpulse=" << tpulse << " pulse2=" << pulse2->Eval(20) << endl;*/

          pulses.push_back(pulse);
          }
   
   hhits->Fill(pulses.size());
   // fill histogram with pulse information

  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);
      
         // loop over all dark pulses
         Double_t amp=0;
         for (j=0; j<pulses.size(); j++) {
	amp += pulses[j]->Eval(t) - pedestal;
	// cout << " jj=" << jj << " npulses=" << pulses.size() << " j=" << j << " amp=" << amp << endl;
                    }
         Double_t fadc = amp + ran_ped;

         hpulse->Fill(amp);
         hpedestal->Fill(ran_ped);
         hsamples->Fill(pedestal+fadc);

         h2_samples->Fill(Vover,pedestal+fadc);
         h2_rate->Fill(dark_rate/1e6,pedestal+fadc);

         // cout << "jj=" << jj << " t=" << t << "ran_ped=" << ran_ped << " fadc=" << fadc << endl;

         if (k == kmax-6 && kk==1) {

         // fill window histogram only with last Vover

             cout << " k=" << k << " 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);
             }
         }

   }     // end over number of events per voltage setting

   }    // end over gain setting.
         
  
   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("Dark hits");
   hpulsew->GetXaxis()->SetTitle("FADC samples");
   hpulsew->GetXaxis()->SetNdivisions(505);
   hpulsew->SetLineColor(2);
   hpulsew->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");
   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();       
  
   TCanvas *c2 = new TCanvas("c2", "c2",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   c2->SetGridy();

   c2->Divide(2,2);
   c2->cd(1);
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);

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

   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);

   hpulse->SetTitle("");
   hpulse->GetXaxis()->SetLabelSize(0.05);
   hpulse->GetXaxis()->SetTitleSize(0.05);
   hpulse->GetYaxis()->SetLabelSize(0.05);
   hpulse->GetYaxis()->SetTitleSize(0.05);
   hpulse->GetYaxis()->SetTitleOffset(1.5);
   hpulse->GetYaxis()->SetTitle("Counts");
   hpulse->GetXaxis()->SetTitle("Pulse only");
   hpulse->GetXaxis()->SetNdivisions(505);
   hpulse->SetLineColor(2);
   hpulse->Draw();    

   c2->cd(3);
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);

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

   c2->cd(4);
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);

   hhits->SetTitle("");
   hhits->GetXaxis()->SetLabelSize(0.05);
   hhits->GetXaxis()->SetTitleSize(0.05);
   hhits->GetYaxis()->SetLabelSize(0.05);
   hhits->GetYaxis()->SetTitleSize(0.05);
   hhits->GetYaxis()->SetTitleOffset(1.5);
   hhits->GetYaxis()->SetTitle("Counts");
   hhits->GetXaxis()->SetTitle("Number of dark hits/window");
   hhits->GetXaxis()->SetNdivisions(505);
   hhits->SetLineColor(2);
   hhits->Draw();   


   TCanvas *c3 = new TCanvas("c3", "c3",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);

   c3->SetGridx();
   c3->SetGridy();

   c3->Divide(2,2);
   c3->cd(1);
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);

   h2_samples->SetTitle("");
   h2_samples->GetXaxis()->SetLabelSize(0.05);
   h2_samples->GetXaxis()->SetTitleSize(0.05);
   h2_samples->GetYaxis()->SetLabelSize(0.05);
   h2_samples->GetYaxis()->SetTitleSize(0.05);
   h2_samples->GetYaxis()->SetTitleOffset(1.5);
   h2_samples->GetYaxis()->SetTitle("Samples");
   h2_samples->GetXaxis()->SetTitle("Vover (V)");
   h2_samples->GetXaxis()->SetNdivisions(505);
   h2_samples->SetLineColor(2);
   h2_samples->Draw("colz");   

    // get mean/sigma from fitting slices

     h2_samples->FitSlicesY();
     TH1D *h2_samples_1 = (TH1D*) gROOT->FindObject("h2_samples_1");
     TH1D *h2_samples_2 = (TH1D*) gROOT->FindObject("h2_samples_2");

     /*for (j=1;j<nbins+1;j++) {
	Float_t mean = h2_samples_1->GetBinContent(j);
	Float_t sigma = h2_samples_2->GetBinContent(j);
	h2_samples_1->SetBinError(j,sigma);
	// cout << " j=" << j << " mean=" << mean << " error=" << sigma << endl;
	}*/

   c3->cd(2);
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);

   Double_t xmin=0;
   Double_t xmax=2.5;
   Double_t ymin=98;
   Double_t ymax=108;

   h2_samples_1->SetTitle(""); 
   // h2_samples_1->SetOptStat(kFALSE);
   h2_samples_1->GetXaxis()->SetRangeUser(xmin,xmax);
   h2_samples_1->GetYaxis()->SetRangeUser(ymin,ymax);
   h2_samples_1->GetXaxis()->SetLabelSize(0.05);
   h2_samples_1->GetXaxis()->SetTitleSize(0.05);
   h2_samples_1->GetYaxis()->SetLabelSize(0.05);
   h2_samples_1->GetYaxis()->SetTitleSize(0.05);
   h2_samples_1->GetYaxis()->SetTitleOffset(1.5);
   h2_samples_1->GetYaxis()->SetTitle("Mean (counts)");
   h2_samples_1->GetXaxis()->SetTitle("Vover (V)");
   h2_samples_1->GetXaxis()->SetNdivisions(505);
   h2_samples_1->SetMarkerStyle(20);
   h2_samples_1->SetMarkerColor(2);
   h2_samples_1->Draw();  

   Double_t xmin=0;
   Double_t xmax=2.5;
   Double_t ymin=1;
   Double_t ymax=3.0;

   c3->cd(3);
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);

   h2_samples_2->SetTitle("");
   // h2_samples_2->SetOptStat(kFALSE);
   h2_samples_2->GetXaxis()->SetRangeUser(xmin,xmax);
   h2_samples_2->GetYaxis()->SetRangeUser(ymin,ymax);
   h2_samples_2->GetXaxis()->SetLabelSize(0.05);
   h2_samples_2->GetXaxis()->SetTitleSize(0.05);
   h2_samples_2->GetYaxis()->SetLabelSize(0.05);
   h2_samples_2->GetYaxis()->SetTitleSize(0.05);
   h2_samples_2->GetYaxis()->SetTitleOffset(1.5);
   h2_samples_2->GetYaxis()->SetTitle("Width (counts)");
   h2_samples_2->GetXaxis()->SetTitle("Vover (V)");
   h2_samples_2->GetXaxis()->SetNdivisions(505);
   h2_samples_2->SetMarkerStyle(20);
   h2_samples_2->SetMarkerColor(2);
   TF1 *fit = new TF1("fit","[0]+[1]*x*x",xmin,xmax);
   h2_samples_2->Fit("fit");  
   h2_samples_2->Draw();    


    c3->cd(4);

   sprintf (string,"Parameter Settings\n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,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.2,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.2,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.2,0.70,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

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

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

   sprintf (string,"Dark rate =%.1f MHz at Vover=%.1f V\n",dark_rate0/1e6,Vover0);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.55,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


   TCanvas *c4 = new TCanvas("c4", "c4",200,10,700,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);

   c4->SetGridx();
   c4->SetGridy();

   c4->Divide(2,2);
   c4->cd(1);
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);

   h2_rate->SetTitle("");
   h2_rate->GetXaxis()->SetLabelSize(0.05);
   h2_rate->GetXaxis()->SetTitleSize(0.05);
   h2_rate->GetYaxis()->SetLabelSize(0.05);
   h2_rate->GetYaxis()->SetTitleSize(0.05);
   h2_rate->GetYaxis()->SetTitleOffset(1.5);
   h2_rate->GetYaxis()->SetTitle("Samples");
   h2_rate->GetXaxis()->SetTitle("Rate (MHz)");
   h2_rate->GetXaxis()->SetNdivisions(505);
   h2_rate->SetLineColor(2);
   h2_rate->Draw("colz");   

    // get mean/sigma from fitting slices

     h2_rate->FitSlicesY();
     TH1D *h2_rate_1 = (TH1D*) gROOT->FindObject("h2_rate_1");
     TH1D *h2_rate_2 = (TH1D*) gROOT->FindObject("h2_rate_2");

     /*for (j=1;j<nbins+1;j++) {
	Float_t mean = h2_rate_1->GetBinContent(j);
	Float_t sigma = h2_rate_2->GetBinContent(j);
	h2_rate_1->SetBinError(j,sigma);
	// cout << " j=" << j << " mean=" << mean << " error=" << sigma << endl;
	}*/

   c4->cd(2);
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);

   Double_t xmin=0;
   Double_t xmax=150;
   Double_t ymin=98;
   Double_t ymax=108;

   h2_rate_1->SetTitle(""); 
   // h2_rate_1->SetOptStat(kFALSE);
   h2_rate_1->GetXaxis()->SetRangeUser(xmin,xmax);
   h2_rate_1->GetYaxis()->SetRangeUser(ymin,ymax);
   h2_rate_1->GetXaxis()->SetLabelSize(0.05);
   h2_rate_1->GetXaxis()->SetTitleSize(0.05);
   h2_rate_1->GetYaxis()->SetLabelSize(0.05);
   h2_rate_1->GetYaxis()->SetTitleSize(0.05);
   h2_rate_1->GetYaxis()->SetTitleOffset(1.5);
   h2_rate_1->GetYaxis()->SetTitle("Mean (counts)");
   h2_rate_1->GetXaxis()->SetTitle("Rate (MHz)");
   h2_rate_1->GetXaxis()->SetNdivisions(505);
   h2_rate_1->SetMarkerStyle(20);
   h2_rate_1->SetMarkerColor(2);
   h2_rate_1->Draw();  

   Double_t xmin=0;
   Double_t xmax=150;
   Double_t ymin=1;
   Double_t ymax=3.0;

   c4->cd(3);
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);

   h2_rate_2->SetTitle("");
   // h2_rate_2->SetOptStat(kFALSE);
   h2_rate_2->GetXaxis()->SetRangeUser(xmin,xmax);
   h2_rate_2->GetYaxis()->SetRangeUser(ymin,ymax);
   h2_rate_2->GetXaxis()->SetLabelSize(0.05);
   h2_rate_2->GetXaxis()->SetTitleSize(0.05);
   h2_rate_2->GetYaxis()->SetLabelSize(0.05);
   h2_rate_2->GetYaxis()->SetTitleSize(0.05);
   h2_rate_2->GetYaxis()->SetTitleOffset(1.5);
   h2_rate_2->GetYaxis()->SetTitle("Width (counts)");
   h2_rate_2->GetXaxis()->SetTitle("Rate (MHz)");
   h2_rate_2->GetXaxis()->SetNdivisions(505);
   h2_rate_2->SetMarkerStyle(20);
   h2_rate_2->SetMarkerColor(2);
   h2_rate_2->Fit(sqroot);
   h2_rate_2->Draw();    



   sprintf (string,"pedestal_width_%02d_%02d.pdf(",layer,T);
   c0->SaveAs(string);
   sprintf (string,"pedestal_width_%02d_%02d.pdf",layer,T);
   c1->SaveAs(string);
   sprintf (string,"pedestal_width_%02d_%02d.pdf",layer,T);
   c2->SaveAs(string);
   sprintf (string,"pedestal_width_%02d_%02d.pdf",layer,T);
   c3->SaveAs(string);
   sprintf (string,"pedestal_width_%02d_%02d.pdf)",layer,T);
   c4->SaveAs(string);

   // save width vs rate histogram.

   sprintf (string,"pedestal_width_c3_%02d_%02d.C",layer,T);
   c3->SaveAs(string);
   sprintf (string,"pedestal_width_c4_%02d_%02d.C",layer,T);
   c4->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;
}
Double_t sqroot_func (Double_t *x, Double_t *par) {

   Double_t parA=par[0];
   Double_t parB=par[1];
   Double_t x1=x[0];

   char string[256];

   Double_t func;

  func = sqrt(parA*parA + parB*parB*x1*x1);
   
   /*sprintf (string,"sqroot_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;
}