void gen_bcal_trig(void)
{
// 
// Generate 4n BCAL slices inputs to the trigger to study thresholds and rates  
//

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

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

   char string[256];
   char filename[80];
   Int_t k,kk,j,jj;
   Int_t nevents=1000;
   
   Double_t mean=2;       // number of dark hits per 100 ns
   Double_t  sigma=sqrt(mean);

   TRandom1 *r = new TRandom1();
   UInt_t  iseed = 23422777519;
   r->SetSeed(iseed);

   Double_t pedestal = 5;
   Double_t pedestal_width = 1.3; 
   // Double_t pedestal = 0;
   // Double_t pedestal_width = 0; 
   Double_t counts_per_pixel=0.3;   // estimate of number of ADC counts per pixel fired
   Int_t nFlash =  1536;     // number of flash channels
   // Int_t nFlash =  15;     // number of flash channels
   Int_t nsigma = 2;     // number of sigma required for selection
   
   // define histograms

   Int_t nbins=25;
   Int_t nbins2=100;

   // compute mean per channel in window

   Double_t tmin=0;
   Double_t tmax=100;

   Int_t sigma_win = 5;
   Int_t mean_win = pedestal+mean*counts_per_pixel/nbins;
   Int_t sig_win = sqrt(pedestal_width*pedestal_width+mean*counts_per_pixel/nbins);
   TH1F *window = new TH1F("window","window distribution",nbins,tmin,tmax);
   TH1F *selection1 = new TH1F("selection1_trig","selection1_trig  distribution",nbins2,mean_win-sigma_win*sig_win,mean_win+sigma_win*sig_win);
   TH1F *selection2 = new TH1F("selection2_trig","selection2_trig  distribution",nbins2,mean_win-sigma_win*sig_win,mean_win+sigma_win*sig_win);

   Int_t sigma_win = 5*sqrt(nFlash);
   Int_t mean_win = nFlash*(pedestal+mean*counts_per_pixel/nbins);
   TH1F *window_trig = new TH1F("window_trig","Summed trigger window distribution",nbins,tmin,tmax);
   TH1F *selection1_trig = new TH1F("selection1_trig","selection1_trig  distribution",nbins2,mean_win-sigma_win*pedestal_width,mean_win+sigma_win*pedestal_width);
   TH1F *selection2_trig = new TH1F("selection2_trig","selection2_trig  distribution",nbins2,mean_win-sigma_win*pedestal_width,mean_win+sigma_win*pedestal_width);
   TH1F *sum_trig = new TH1F("sum_trig","sum_trig  distribution",nbins2,nbins*(mean_win-sigma_win*pedestal_width/sqrt(nbins)),nbins*(mean_win+sigma_win*pedestal_width/sqrt(nbins)));

   // define normalized pulse shape
  // assume pulse is located at x=0 (just below the peak)

   Double_t xmin=-20;
   Double_t xmax=120;
   // parameters obtained from fit_bcal_pulse and trial/error
   // Double_t pulse_mu=0.04;  // units of 1/x   eyeball fit - poor recovery to baseline
   // Double_t pulse_sigma=3.5;     // units of x eyeball fit - matches fast rise time
   Double_t pulse_mu=0.058;  // units of 1/x   actual fit  - good recovery to baseline
   Double_t pulse_sigma=5.91;     // units of x actualfit  - slower rise time
   Double_t gain=counts_per_pixel;   

   TF1 *pulse = new TF1("pulse_func",pulse_func,xmin,xmax,3);
   pulse->SetParameters(pulse_mu,pulse_sigma,gain);
   Double_t pulse_sum = pulse->Integral(xmin,xmax);
   // printf ("pulse_sum=%f \n",pulse_sum);

#if 0
   //
   TCanvas *c10 = new TCanvas("c10","c10 gen_bcal_trig",200,10,700,700);
   c10->SetBorderMode(0);
   c10->SetFillColor(0);
   // c10->SetGridx();
   // c10->SetGridy();
   // c10->SetLogy();
   pulse->Draw("");
#endif

   Double_t sum=0;
   Double_t delta=4;

  // loop over events

  for (kk=0;kk<nevents;kk++) {
 
  // loop over number of flash ADCs
  // first reset window and continue loop

  window_trig->Reset();
  selection1_trig->Reset();
  selection2_trig->Reset();
  selection1->Reset();
  selection2->Reset();

  for (k=0;k<nFlash;k++) {

    // first reset window and continue loop

   window->Reset();

  // generate window distribution. First generate pedestal

  for (j=0;j<nbins;j++) {
         Double_t x=tmin + (tmax-tmin)*(j+0.5)/nbins;
         Double_t y= r->Gaus(pedestal,pedestal_width);
         window->Fill(x,y);
         }


  // Dark pulse generation: need to generate x beyond the range of the histogram to properly fill the phase space with pulses

  Double_t delta_win = 50;
  Double_t tmin_win = tmin - delta_win;
  Double_t tmax_win = tmax + delta_win;

   // Generate number of dark pulses
   Int_t ndark = r->Poisson(mean*(tmax_win-tmin_win)/100);

   // printf ("\nGenerate ndark=%d hits\n",ndark);

   for (j=0;j<ndark;j++) {

       Double_t x = tmin_win + (tmax_win - tmin_win)*r->Uniform();
      
       // fill one pulse per pixel fired at location x

       for (jj=0;jj<25;jj++) {
          Double_t xpulse = jj*delta + xmin;
          if (x + xpulse < tmin | | x + xpulse > tmax) continue;
          Double_t y = delta*pulse->Eval(xpulse);
          sum += y;
          // printf ("x=%f xpulse=%f y=%f sum=%f \n",x,xpulse,y,sum);
          window->Fill(x+xpulse,y);
          }


       }


    // input data into selection histograms from single FADC window histogram


   // Create histogram of entries (projection onto y axis)

      Int_t ndim1 = window->GetNbinsX();
      // printf ("\nhist_read: ndim1=%d, xlo=%f, width=%f\n\n",ndim1,xlo,width);
     Double_t average=0;
     Double_t rms2=0;

      for(j=0;j<ndim1;j++) {
          Double_t xbin = window->GetBinCenter(j+1);
          Double_t content = window->GetBinContent(j+1);
          selection1->Fill(content);
          average += content;
          rms2 += content*content;
          // printf ("selection1_trig, j=%d, xbin=%f, content=%f, average=%f\n",j,xbin,content,average);
         }
     average /= ndim1;
     rms2 /= ndim1;
     Double_t rms = sqrt(rms2 - average*average);
     // printf ("pedestal trig=%f, pedestal_width=%f, average=%f, rms=%f\n",pedestal,pedestal_width,average,rms);

     // get selection 2 based on 

      for(j=0;j<ndim1;j++) {
          Double_t xbin = window->GetBinCenter(j+1);
          Double_t content = window->GetBinContent(j+1);
          if (content > average + nsigma*rms) {
          // if (content > pedestal + nsigma*pedestal_width)) {
                  selection2->Fill(content);
                  // printf ("selection2_trig, j=%d, xbin=%f, content=%f\n",j,xbin,content);
                  }
          }


     // sum windows from individual FADCs into a total trigger sum

    *window_trig = *window_trig +  *window;
    }

   // Create histogram of entries (projection onto y axis)

      Int_t ndim1 = window_trig->GetNbinsX();
      // printf ("\nhist_read: ndim1=%d, xlo=%f, width=%f\n\n",ndim1,xlo,width);
     Double_t average=0;
     Double_t rms2=0;

      for(j=0;j<ndim1;j++) {
          Double_t xbin = window_trig->GetBinCenter(j+1);
          Double_t content = window_trig->GetBinContent(j+1);
          selection1_trig->Fill(content);
          average += content;
          rms2 += content*content;
          // printf ("selection1_trig, j=%d, xbin=%f, content=%f, average=%f\n",j,xbin,content,average);
         }
     average /= ndim1;
     rms2 /= ndim1;
     Double_t rms = sqrt(rms2 - average*average);
     // printf ("pedestal trig=%f, pedestal_width=%f, average=%f, rms=%f\n",pedestal,pedestal_width,average,rms);

     // get selection 2 based on 

      for(j=0;j<ndim1;j++) {
          Double_t xbin = window_trig->GetBinCenter(j+1);
          Double_t content = window_trig->GetBinContent(j+1);
          if (content > average + nsigma*rms) {
          // if (content > pedestal + nsigma*pedestal_width)) {
                  selection2_trig->Fill(content);
                  // printf ("selection2_trig, j=%d, xbin=%f, content=%f\n",j,xbin,content);
                  }
          }


   // fill histogram with integral contents over events

  sum_trig->Fill(window_trig->Integral());
      
  }

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

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

   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);
   window->SetTitle("");
   window->GetYaxis()->SetRangeUser(0,2*window->GetMaximum());   
   window->GetXaxis()->SetTitleSize(0.04);
   window->GetYaxis()->SetTitleSize(0.04);
   window->GetYaxis()->SetTitleOffset(1.5);
   window->GetXaxis()->SetTitle("Time (4n bins)");
   window->GetYaxis()->SetTitle("One Flash (FADC Counts)");
   window->GetXaxis()->SetNdivisions(505);
   window->Draw();
  
   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   selection1->SetTitle("");
   // selection1->GetYaxis()->SetRangeUser(0,2*pedestal);   
   selection1->GetXaxis()->SetTitleSize(0.04);
   selection1->GetYaxis()->SetTitleSize(0.04);
   selection1->GetYaxis()->SetTitleOffset(1.5);
   selection1->GetXaxis()->SetTitle("Amplitude (mV)");
   selection1->GetYaxis()->SetTitle("Entries");
   selection1->GetXaxis()->SetNdivisions(505);
   selection1->Draw();
  
   sprintf (string,"Number of Flashes= %d\n",nFlash);
   // printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

  
   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   selection2->SetTitle("");
   // selection2->GetYaxis()->SetRangeUser(0,2*pedestal);   
   selection2->GetXaxis()->SetTitleSize(0.04);
   selection2->GetYaxis()->SetTitleSize(0.04);
   selection2->GetYaxis()->SetTitleOffset(1.5);
   selection2->GetXaxis()->SetTitle("Amplitude (mV)");
   selection2->GetYaxis()->SetTitle("Entries");
   selection2->GetXaxis()->SetNdivisions(505);
   selection2->Draw();
  
   sprintf (string,"Selection 2 based on %d#sigma\n",nsigma);
   // printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();



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

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

   c2->cd(1);
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   window_trig->SetTitle("");
   window_trig->GetYaxis()->SetRangeUser(0,2*window_trig->GetMaximum());   
   window_trig->GetXaxis()->SetTitleSize(0.04);
   window_trig->GetYaxis()->SetTitleSize(0.04);
   window_trig->GetYaxis()->SetTitleOffset(1.5);
   window_trig->GetXaxis()->SetTitle("Time (4n bins)");
   window_trig->GetYaxis()->SetTitle("Sum All Flashes (FADC Counts)");
   window_trig->GetXaxis()->SetNdivisions(505);
   window_trig->Draw();
  
   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   selection1_trig->SetTitle("");
   // selection1_trig->GetYaxis()->SetRangeUser(0,2*pedestal);   
   selection1_trig->GetXaxis()->SetTitleSize(0.04);
   selection1_trig->GetYaxis()->SetTitleSize(0.04);
   selection1_trig->GetYaxis()->SetTitleOffset(1.5);
   selection1_trig->GetXaxis()->SetTitle("Amplitude (mV)");
   selection1_trig->GetYaxis()->SetTitle("Entries");
   selection1_trig->GetXaxis()->SetNdivisions(505);
   selection1_trig->Draw();
  
   sprintf (string,"Number of Flashes= %d\n",nFlash);
   // printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


  
   c2->cd(3);
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   selection2_trig->SetTitle("");
   // selection2_trig->GetYaxis()->SetRangeUser(0,2*pedestal);   
   selection2_trig->GetXaxis()->SetTitleSize(0.04);
   selection2_trig->GetYaxis()->SetTitleSize(0.04);
   selection2_trig->GetYaxis()->SetTitleOffset(1.5);
   selection2_trig->GetXaxis()->SetTitle("Amplitude (mV)");
   selection2_trig->GetYaxis()->SetTitle("Entries");
   selection2_trig->GetXaxis()->SetNdivisions(505);
   selection2_trig->Draw();
  
   sprintf (string,"Selection 2 based on %d#sigma\n",nsigma);
   // printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();

   c2->cd(4);
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   sum_trig->SetTitle("");
   // sum_trig->GetYaxis()->SetRangeUser(0,2*pedestal);   
   sum_trig->GetXaxis()->SetTitleSize(0.04);
   sum_trig->GetYaxis()->SetTitleSize(0.04);
   sum_trig->GetYaxis()->SetTitleOffset(1.5);
   sum_trig->GetXaxis()->SetTitle("Amplitude (mV)");
   sum_trig->GetYaxis()->SetTitle("Sum Over Window all Flashes(FADC Counts)");
   sum_trig->GetXaxis()->SetNdivisions(505);
   sum_trig->Draw();
  
   sprintf (string,"Sum all based on %d events\n",nevents);
   // printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.04);
   t1->Draw();


   sprintf(filename,"gen_bcal_trig_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"gen_bcal_trig_c1.png");
   c1->SaveAs(filename);
   sprintf(filename,"gen_bcal_trig_c2.pdf");
   c2->SaveAs(filename);
   sprintf(filename,"gen_bcal_trig_c2.png");
   c2->SaveAs(filename);

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

   // pulse function with Gaussian convoluted with exponential
  // the distribution is normalized to one.

   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 arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);
   
   amplitude = gain * amplitude;
   
   /*sprintf (string,"amplitude=%f mu=%f sigma=%f gain=%f arg=%f\n",amplitude,mu,sigma,gain,arg);
   printf ("string=%s",string);*/

   return amplitude;
}