void apd_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(kFALSE);
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;

   // 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->Divide(2,2);
   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->SetLogy();

   c1->cd(1);
   c1_1->SetGridx();
   c1_1->SetGridy();
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);
   c1_1->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=250;
   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.90;
   Double_t coverage=0.90;
   Double_t egain0=250;
   Double_t egain=egain0;
   Double_t atten=exp(-length/(2*lambda));
   Double_t gain=ng*lg*pde*coverage*sgain*egain;
   Double_t gain0=ng*lg*pde*coverage*sgain*egain;
   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

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.35,0.86,0.85,0.9);
   dummy->SetMarkerColor(1);
   dummy->SetMarkerStyle(21);
   t1 = new TLatex(0.15,0.91,"Apd: Min pulse of max channel");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw(); 
   dummy->SetTitle("");
   dummy->GetXaxis()->SetRangeUser(xmin,xmax);
   dummy->GetYaxis()->SetRangeUser(ymin,ymax);
   dummy->GetXaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleSize(0.04);
   dummy->GetYaxis()->SetTitleOffset(1.5);
   dummy->GetXaxis()->SetTitle("Time (ns)");
   dummy->GetYaxis()->SetTitle("Amplitude (mV)");
   dummy->Draw("Ap");
   dummy->GetXaxis()->SetNdivisions(505);

   // plot pulse function

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


   Double_t f1=2;    
   pulse->SetParameter(2,f1*Eg*gain*ratio90to15/atten);
   printf ("amplification=%e, gain=%e\n",f1*Eg*gain*ratio90to15/atten,gain);
   Double_t mean_npe = f1*Eg*ng_per_GeV*lg*pde*coverage/atten;
   pulse->SetParameter(3,mean_npe);

   // units are 50 Ohm mV x ns
   Double_t int1 = pulse->Integral(xmin,xmax)*1e-3*1e-9/termination;
   sprintf (string,"#mu=%.3f, #sigma=%.2f, Q=%.2e C\n",mu,sigma,int1);
   // sprintf(string,"Q: %.2e pC\n",int1*1e12);
   // printf ("int1a %s",string);
   TF1 *pulse1 = pulse->DrawCopy("sameC");
   leg->AddEntry(pulse1,string,"l");
   leg->SetLineColor(1);
   leg->SetTextSize(0.025);
   leg->Draw();

   Double_t f2=0.5;
   pulse->SetParameter(2,f2*gain);
   Double_t mean_npe = f2*Eg*ng_per_GeV*lg*pde*coverage;
   pulse->SetParameter(3,mean_npe);
   pulse->SetLineColor(2);
   TF1 *pulse2 = pulse->DrawCopy("sameC");
   Double_t int2 = pulse2->Integral(xmin,xmax)*1e-3*1e-9/termination;

   Double_t f3=0.04;

   // allow for fluctuations at small number of photoelectrons

   TRandom1 *r = new TRandom1();
   r->SetSeed(0);
   Double_t mean_npe = f3*Eg*atten*ng_per_GeV*lg*pde*coverage;
   Double_t npe = (Double_t) r->Poisson(mean_npe);

   TF1 *pulse3 = 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   
   pulse3->SetParameter(0,mu);    
   pulse3->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);
      }
   }

   // adjust gain for actual number of photoelecrons

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

   pulse3->SetParameter(2,f3*Eg*gain*atten);
   pulse3->SetParameter(3,npe);
   for (j=0;j<npe && j<maxnpe;j++) {
      pulse3->SetParameter(4+j,texp[j]);
      printf ("j=%d, texp[j]=%f\n",j, texp[j]);
   }

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

   sprintf(string,"Sensor gain          = %.2e\n",sgain0);
   t1 = new TLatex(0.50,0.79,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf(string,"E_{#gamma}                           = %.2e GeV\n",Eg);
   t1 = new TLatex(0.50,0.76,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw(); 
   sprintf(string,"n_{#gamma}/E_{#gamma}                       = %.2e\n",ng_per_GeV);
   t1 = new TLatex(0.50,0.73,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Light Collection    = %.2e\n",lg);
   t1 = new TLatex(0.50,0.70,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Sensor Coverage  = %.2e\n",coverage);
   t1 = new TLatex(0.50,0.67,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"PDE                        = %.2e\n",pde);
   t1 = new TLatex(0.50,0.64,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Electronic gain     = %.2e\n",egain);
   t1 = new TLatex(0.50,0.61,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Amplif Factor       = %.2e\n",gain0);
   t1 = new TLatex(0.50,0.58,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Atten  Factor        = %.2e\n",atten);
   t1 = new TLatex(0.50,0.55,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"15deg   Factor        = %.2e\n",ratio90to15);
   t1 = new TLatex(0.50,0.52,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();

   Double_t th=15;
   TLine *thresh = new TLine(xmin,th,xmax,th);
   thresh->SetLineColor(2);
   thresh->Draw();

   sprintf(string,"E_{#gamma}=%.1f GeV 15deg / atten\n",f1*Eg);
   t1 = new TLatex(0.5,0.45,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"E_{#gamma}=%.1f GeV\n",f2*Eg);
   t1 = new TLatex(0.5,0.42,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"E_{#gamma}=%.3f GeV * atten\n",f3*Eg);
   t1 = new TLatex(0.5,0.39,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Thresh=%.0f mV\n",th);
   t1 = new TLatex(0.65,0.2,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();

   // Create FADC samples from each pulse.
   
   c1->cd(2);
   c1_2->SetGridx();
   c1_2->SetGridy();
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   // c1_2->SetLogy();

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

   if (nsamples > msamples) {
      printf ("*** Apd_pulse - nsamples=%d, msamples=%d\n",nsamples,msamples);
      nsamples = msamples;
      }
   Double_t adc1[maxsamples];
   Double_t adc2[maxsamples];
   Double_t adc3[maxsamples];
   Int_t fadc1[maxsamples];
   Int_t fadc2[maxsamples];
   Int_t fadc3[maxsamples];
   Int_t tbins[maxsamples];
   Double_t sum1=0;
   Double_t sum2=0;
   Double_t sum3=0;

   for (j=0;j<nsamples;j++) {
      Double_t t = xmin + j*tsample;
      tbins[j] = j;
      adc1[j] = pulse1->Eval(t);
      adc2[j] = pulse2->Eval(t);
      adc3[j] = pulse3->Eval(t);
      fadc1[j] = maxbits*adc1[j]/Vmax;
      sum1 = sum1 + fadc1[j];
      fadc2[j] = maxbits*adc2[j]/Vmax;
      sum2 = sum2 + fadc2[j];
      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,adc1[j],fadc1[j]);
      // printf ("string=%s",string);
      }

   TGraph *p1 = new TGraph (nsamples,tbins,fadc1);
   t1 = new TLatex(0.20,0.91,"Apd: Pulse max channel");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   // t1->Draw(); 
   p1->SetTitle("");
   p1->GetXaxis()->SetRangeUser(0,nsamples);
   p1->GetYaxis()->SetRangeUser(0,6000*egain/egain0);
   p1->GetXaxis()->SetTitleSize(0.04);
   p1->GetYaxis()->SetTitleSize(0.04);
   p1->GetYaxis()->SetTitleOffset(2);
   p1->GetXaxis()->SetTitle("Time stamp (4 ns)");
   p1->GetYaxis()->SetTitle("FADC (counts)");
   p1->SetMarkerColor(1);
   p1->SetMarkerStyle(21);
   p1->SetMarkerSize(0.5);
   p1->GetXaxis()->SetNdivisions(505);
   p1->Draw("Ap");

   sprintf(string,"E_{#gamma}=%.1f GeV *15deg /atten, sum=%.0f\n",f1*Eg,sum1);
   t1 = new TLatex(0.3,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"FADC: 0-%.0f mV, 0-%d\n",Vmax,maxbits);
   t1 = new TLatex(0.3,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"Q: %.2e pC\n",int1*1e12);
   printf ("int1 %s",string);
   t1 = new TLatex(0.3,0.70,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   
   
   c1->cd(3);
   c1_3->SetGridx();
   c1_3->SetGridy();
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   // c1_3->SetLogy();


   TGraph *p2 = new TGraph (nsamples,tbins,fadc2);
   t1 = new TLatex(0.20,0.91,"Apd: Pulse max channel");
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   // t1->Draw(); 
   p2->SetTitle("");
   p2->GetXaxis()->SetRangeUser(0,nsamples);
   p2->GetYaxis()->SetRangeUser(0,200*egain/egain0);
   p2->GetXaxis()->SetTitleSize(0.04);
   p2->GetYaxis()->SetTitleSize(0.04);
   p2->GetYaxis()->SetTitleOffset(2);
   p2->GetXaxis()->SetTitle("Time stamp (4 ns)");
   p2->GetYaxis()->SetTitle("FADC (counts)");
   p2->SetMarkerColor(2);
   p2->SetMarkerStyle(21);
   p2->SetMarkerSize(0.5);
   p2->GetXaxis()->SetNdivisions(505);
   p2->Draw("Ap");
   sprintf(string,"E_{#gamma}=%.1f GeV, sum=%.0f\n",f2*Eg,sum2);
   t1 = new TLatex(0.3,0.8,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"FADC: 0-%.0f mV, 0-%d\n",Vmax,maxbits);
   t1 = new TLatex(0.3,0.75,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"Q: %.2e pC\n",int2*1e12);
   printf ("int2 %s",string);
   t1 = new TLatex(0.3,0.70,string);
   t1->SetTextColor(2);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   
   c1->cd(4);
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);
   // c1_4->SetLogy();


   TGraph *p3 = new TGraph (nsamples,tbins,fadc3);
   t1 = new TLatex(0.20,0.91,"Apd: 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,10*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.3,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.3,0.75,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

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

   sprintf(string,"(#approx %.1f mip equivalent)\n",frac_mip);
   t1 = new TLatex(0.3,0.70,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();
   sprintf(string,"#LTn_{pe}#GT=%.1f, npe=%d\n",mean_npe,npe);
   t1 = new TLatex(0.3,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.3,0.60,string);
   t1->SetTextColor(4);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


   Int_t iEg = 1000*(Eg+0.0005);
   sprintf(filename,"apd_pulse_%dMeV.eps",iEg);
   c1->SaveAs(filename);
   sprintf(filename,"apd_pulse_%dMeV.gif",iEg);
   c1->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 arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);
   

   // 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 = par[4+j] + timeoffset;
         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;
         }
      }

   // convert to voltage scale (mV).

   amplitude = gain * amplitude * conversion;
   return amplitude;
}