void small_pulse(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//

#include <TMath.h>
#include <TRandom.h>

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

   char string[256];
   char filename[80];
   Int_t j,jj;
   
   #define npts 3;
   #define maxsamples 100;
   #define maxnpe 10;
   #define npulses 1000;

   // input data set first point (actually -100) to 100 so that it is not plotted.

   // use offset on x-axis to make data visible and give negative y-values for reference (not plotted)

   Double_t dummyx[npts]={-100,0,160};
   Double_t dummyy[npts]={-1,-1,-1};

   //
   TCanvas *c1 = new TCanvas("c1","c1 Sensor pulses",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetGridx();
   c1->SetGridy();
   // c1->SetLogy();


   // use energy deposited in cell with max energy.

   Double_t lg0=0.8;
   Double_t length=390;
   Double_t lambda=280;
   Double_t sgain0=1e6;
   Double_t sgain=sgain0;
   Double_t Eg=1;
   Double_t ng_per_GeV=1814;
   Double_t ng=ng_per_GeV*Eg;
   Double_t lg=0.8/lg0;     // relative to lg in beamtest
   Double_t pde=0.2;
   // Double_t pde=0.4;
   Double_t coverage=0.77;
   Double_t egain0=10;
   // Double_t egain0=5;
   Double_t egain=egain0*2/6.3;
   Double_t atten=exp(-length/(2*lambda));
   Double_t gain0=ng*lg*pde*coverage*sgain0*egain;
   Double_t gain=gain0;
   Double_t ratio90to15=3;
   Double_t nphot_mip=44*lg*coverage/0.21;

   Double_t xmin=-20;
   Double_t xmax=120;
   Double_t ymin=gain*1e-8/4;
   // Double_t ymin=0;
   Double_t ymax=gain*1e-5;
   Double_t termination=50;  // 50 Ohms


   // Create FADC samples from each pulse.

   // Double_t f3=0.040;
   Double_t f3=0.060;
   Double_t mean_npe = f3*Eg*atten*ng_per_GeV*lg*pde*coverage;
   Double_t nphotons = f3*Eg*atten*ng_per_GeV*lg*coverage;
   Double_t frac_mip = nphotons/nphot_mip;
   printf("nphotons=%f, mean_npe=%f\n",nphotons,mean_npe);

   // loop over various pulse shapes

   TH2F *digisum = new TH2F ("digisum","Digital sum vs npe",20,0,20,100,0,300);
   TProfile *profsum = new TProfile ("profsum","Profile sum vs npe",20,0,20,"s");

   for (jj=0;jj<npulses;jj++) {

   TRandom1 *r = new TRandom1();
   // UInt_t seed=1111111111;
   r->SetSeed(0);
   Double_t npe = (Double_t) r->Poisson(mean_npe);

   // plot pulse function

   TF1 *pulse = new TF1("pulse_func",pulse_func,xmin,xmax,4+npe);
   Double_t mu=1/16.;    // unis of 1/x
   // Double_t sigma=8;     // units of x 
   Double_t sigma=3;     // units of x   
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);


   TRandom1 *rexp = new TRandom1();
   rexp->SetSeed(0);
   Double_t texp[maxnpe];

   if (npe<maxnpe) {

      for (j=0;j<npe;j++) {
          texp[j] = rexp->Exp(1/mu);

         // include max on exponential list

         /*while (texp[j] > 4*1/mu) {
            Double_t texp[j] = rand->Exp(1/mu);
         }*/
      }
   }

   // adjust gain for actual number of photoelecrons

   sgain = sgain0*npe/mean_npe;
   gain = gain0*npe/mean_npe;

   pulse->SetParameter(2,f3*gain*atten/Eg);
   // printf ("pde=%f, egain=%f, f3*gain*atten/Eg=%f\n",pde,egain,f3*gain*atten/Eg);

   pulse->SetParameter(3,npe);
   for (j=0;j<npe && j<maxnpe;j++) {
      pulse->SetParameter(4+j,texp[j]);
      // printf ("j=%d, texp[j]=%f\n",j, texp[j]);
   }

   pulse->SetLineColor(4);
   Double_t int3 = pulse->Integral(xmin,xmax)*1e-3*1e-9/termination;
   pulse->Draw("");

   Double_t tsample=4;
   Int_t nsamples = (xmax-xmin)/tsample;
   Double_t Vmax=2000;    // 2000 mV.
   Int_t msamples=maxsamples;
   Double_t maxbits=4095;

   if (nsamples > msamples) {
      printf ("*** small_pulse - nsamples=%d, msamples=%d\n",nsamples,msamples);
      nsamples = msamples;
      }
   Double_t adc3[maxsamples];
   Int_t fadc3[maxsamples];
   Int_t tbins[maxsamples];
   Double_t sum3=0;

   for (j=0;j<nsamples;j++) {
      Double_t t = xmin + j*tsample;
      tbins[j] = j;
      adc3[j] = pulse->Eval(t);
      fadc3[j] = maxbits*adc3[j]/Vmax;
      sum3 = sum3 + fadc3[j];
      sprintf (string,"nsamples=%d, j=%d, t=%f, adc=%f, fadc=%d\n",nsamples,j,t,adc3[j],fadc3[j]);
      // printf ("string=%s",string);
      }

   // printf ("jj=%d, npe=%f\n",jj,npe);
   digisum->Fill(npe,sum3);
   profsum->Fill(npe,sum3);
   }

   //
   TCanvas *c2 = new TCanvas("c2","c2 Sensor pulses",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->SetGridx();
   c2->SetGridy();
   // c2->SetLogy();

   TGraph *p3 = new TGraph (nsamples,tbins,fadc3);
   t1 = new TLatex(0.20,0.91,"SensL: Pulse max channel");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   // t1->Draw(); 
   p3->SetTitle("");
   p3->GetXaxis()->SetRangeUser(0,nsamples);
   // p3->GetYaxis()->SetRangeUser(0,50*egain/egain0);
   p3->GetXaxis()->SetTitleSize(0.04);
   p3->GetYaxis()->SetTitleSize(0.04);
   p3->GetYaxis()->SetTitleOffset(2);
   p3->GetXaxis()->SetTitle("Time stamp (4 ns)");
   p3->GetYaxis()->SetTitle("FADC (counts)");
   p3->SetMarkerColor(4);
   p3->SetMarkerStyle(21);
   p3->SetMarkerSize(0.5);
   p3->GetXaxis()->SetNdivisions(505);
   p3->Draw("Ap");
   sprintf(string,"E_{#gamma}=%.3f GeV *atten, sum=%.0f\n",f3*Eg,sum3);
   t1 = new TLatex(0.45,0.8,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"FADC: 0-%.0f mV, 0-%d\n",Vmax,maxbits);
   t1 = new TLatex(0.45,0.75,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


   sprintf(string,"(#approx %.1f mip equivalent)\n",frac_mip);
   t1 = new TLatex(0.45,0.70,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"n_{pe}=%.1f\n",npe);
   t1 = new TLatex(0.45,0.65,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"Q: %.2e pC\n",int3*1e12);
   printf ("int3 %s",string);
   t1 = new TLatex(0.45,0.60,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   
   //
   TCanvas *c3 = new TCanvas("c3","c3 Digital variation in pulse sum",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   c3->SetGridx();
   c3->SetGridy();
   // c3->SetLogy();
   digisum->Draw("box");

   //
   TCanvas *c4 = new TCanvas("c4","c4 Digital variation in pulse sum",200,10,700,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   c4->SetGridx();
   c4->SetGridy();
   // c4->SetLogy();
   profsum->Draw();

   Int_t ipe = 100*(mean_npe);
   // sprintf(filename,"small_pulse_%dMeV.eps",ipe);
   sprintf(filename,"small_pulse_c1_%d_%.0f.eps",ipe,sigma);
   c1->SaveAs(filename);
   sprintf(filename,"small_pulse_c1_%d_%.0f.gif",ipe,sigma);
   c1->SaveAs(filename);
   sprintf(filename,"small_pulse_c2_%d_%.0f.eps",ipe,sigma);
   c2->SaveAs(filename);
   sprintf(filename,"small_pulse_c2_%d_%.0f.gif",ipe,sigma);
   c2->SaveAs(filename);
   sprintf(filename,"small_pulse_c3_%d_%.0f.eps",ipe,sigma);
   c3->SaveAs(filename);
   sprintf(filename,"small_pulse_c3_%d_%.0f.gif",ipe,sigma);
   c3->SaveAs(filename);
   sprintf(filename,"small_pulse_c4_%d_%.0f.eps",ipe,sigma);
   c4->SaveAs(filename);
   sprintf(filename,"small_pulse_c4_%d_%.0f.gif",ipe,sigma);
   c4->SaveAs(filename);

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

   Double_t timeoffset=10;
   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t gain=par[2];
   Double_t npe=par[3];
   Double_t x1=x[0]-timeoffset;
   Double_t pi=3.14159;

   Int_t j;
   char string[256];
   Double_t conversion=8e-6;   // use mV as default

   #define maxnpe 10;
   Double_t time[maxnpe];

   if (npe<maxnpe) {
      for (j=0;j<npe;j++) {
         time[j] = par[4+j] + timeoffset;
         }
   }

   Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);

   // convert to voltage scale (mV).
   
   amplitude = gain * amplitude * conversion;

   // sprintf (string,"amplitude=%f mu=%f sigma=%f arg=%f, npe=%f\n",amplitude,mu,sigma,arg,npe);
   // printf ("string=%s",string);

   // consider fluctuations

   if (npe < maxnpe) {
      amplitude = 0;

      for (j=0;j<npe;j++) {
      
         // compute function as sum of Gaussians
         Double_t t = time[j];
         Double_t amp1 = gain*conversion*exp(-(t-x1)*(t-x1)/(2*sigma*sigma))/sqrt(2*pi*sigma*sigma);
         amplitude = amp1 + amplitude;

         // printf ("x1=%f, j=%d, t=%f, sigma=%f, 1/mu=%f, amp1=%f, amplitude=%f\n",x1,j,t,sigma,1/mu,amp1,amplitude);
         }

      // normalize function 

      if (npe > 0) {
         amplitude = amplitude/npe;
         }
      else {
         amplitude = 0;
         }
      }

   return amplitude;
}