void gen_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;
   Int_t ngen=10000;

   // 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->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetLogy();

   // use energy deposited in cell with max energy.

   Double_t sgain=1e6;
   Double_t Eg=1;
   Double_t ng_per_GeV=1814;
   Double_t ng=ng_per_GeV*Eg;
   Double_t lg=0.8;
   Double_t pde=0.2;
   Double_t coverage=0.77;
   Double_t egain=10;
   Double_t gain=ng*lg*pde*coverage*sgain*egain;

   Double_t xmin=-20;
   Double_t xmax=120;
   Double_t ymin=gain*1e-8;
   Double_t ymax=gain*1e-6; // log
   // Double_t ymax=gain*1e-6/2; lin
   Double_t termination=50;  // 50 Ohms

   // dummy to draw axes

   TGraph *dummy = new TGraph (npts,dummyx,dummyy);
   TLegend *leg = new TLegend(0.4,0.7,0.85,0.9);
   dummy->SetMarkerColor(1);
   dummy->SetMarkerStyle(21);
   t1 = new TLatex(0.20,0.91,"Shapes");
   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,3);
   Double_t mu=1/8.;    // unis of 1/x
   Double_t sigma=5;     // units of x    
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);    
   pulse->SetParameter(2,gain);

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

   Double_t mu=1./16.;    
   Double_t sigma=5;    
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);   
   pulse->SetParameter(2,gain);
   pulse->SetLineColor(2);
   // units are 50 Ohm mV x ns
   Double_t sum = pulse->Integral(xmin,xmax)*1e-3*1e-9/termination;
   sprintf (string,"#mu=%.3f, #sigma=%.2f, Q=%.2e C\n",mu,sigma,sum);
   TF1 *pulse2 = pulse->DrawCopy("sameC");
   leg->AddEntry(pulse2,string,"l");

   Double_t mu=1/32.;    
   Double_t sigma=5;    
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);   
   pulse->SetParameter(2,gain);
   pulse->SetLineColor(4);
   // units are 50 Ohm mV x ns
   Double_t sum = pulse->Integral(xmin,xmax)*1e-3*1e-9/termination;
   sprintf (string,"#mu=%.3f, #sigma=%.2f, Q=%.2e C\n",mu,sigma,sum);
   TF1 *pulse3 = pulse->DrawCopy("sameC");
   leg->AddEntry(pulse3,string,"l");

   Double_t mu=1/64.;    
   Double_t sigma=5;    
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);   
   pulse->SetParameter(2,gain);
   pulse->SetLineColor(3);
   // units are 50 Ohm mV x ns
   Double_t sum = pulse->Integral(xmin,xmax)*1e-3*1e-9/termination;
   sprintf (string,"#mu=%.3f, #sigma=%.2f, Q=%.2e C\n",mu,sigma,sum);
   TF1 *pulse4 = pulse->DrawCopy("sameC");
   leg->AddEntry(pulse4,string,"l");
   leg->Draw();

   sprintf(string,"Sensor gain          = %.2e\n",sgain);
   t1 = new TLatex(0.50,0.65,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.62,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.59,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.56,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.53,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"PDE                        = %.2e\n",pde);
   t1 = new TLatex(0.50,0.50,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.47,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
   sprintf(string,"Overall Amplif       = %.2e\n",gain);
   t1 = new TLatex(0.50,0.44,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.025);
   t1->Draw();
  

   Int_t iEg = 1000*(Eg+0.0005);
   sprintf(filename,"gen_pulse_%dMeV_log.eps",iEg);
   c1->SaveAs(filename);
   sprintf(filename,"gen_pulse_%dMeV_log.gif",iEg);
   c1->SaveAs(filename);

}
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 x1=x[0];

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

   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\n",amplitude,mu,sigma,arg);
   printf ("string=%s",string);*/

   return amplitude;
}