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 {
        
        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,"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;
}
void plot_waveform_integral(void)
{

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
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;
    
    
    ofstream listfile;
    listfile.open("plot_waveform_integral.list");
    
    

   // histogram parameters

  Int_t nbins=100;
   Double_t mu=0.0787;    // unis of 1/x
   Double_t sigma=0.865;     // units of x 
   Double_t gain=350;   
   Double_t t0=39.5;
   Double_t pedestal=2;
    
   // define histogram to hold pulses

  Double_t xmin=0;
  Double_t xmax=100;
    
  TH1F *h_wave= new TH1F("h_wave","Digitized waveform", nbins,xmin,xmax);

   // define pulse function: time units are FADC samples.

   Double_t xpmin=0;
   Double_t xpmax=200;

   TF1 *pulse = new TF1("pulse_func",pulse_func,xpmin,xpmax,5);
   pulse->SetParameters(mu,sigma,gain,t0,pedestal);
   pulse->SetParNames("mu","sigma","gain","t0","pedestal");

   // fill histogram with pulse information

  Double_t Rup = 0.5;     // up reflection coefficient
  Double_t Rdn = 0.5;     // up reflection coefficient
  Double_t L = 390;        // length of BCAL module
  Double_t veff = 17;      // effective velocity in module
  Double_t xpos = -0*L;   // position of interaction
  Double_t lambda = 436;   // GlueX-doc-2049

  for (jj=0;jj<nbins;jj++) {
      Double_t time0 = xmin + (xmax-xmin)*(jj+0.5)/nbins;
      Double_t time = -xpos/veff + time0;
      // Double_t atten = exp(-(L/2+xpos)/lambda);
      Double_t atten = 1;
      // Double_t x = time/4 + t0;
      Double_t x = time + t0 - 20;
      Double_t y = atten*(pulse->Eval(x) - pedestal);
      h_wave->Fill(time0,y);
      
         }
         
  
   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);
   gPad->SetBorderMode(0);
   gPad->SetFillColor(0);

   h_wave->SetTitle("");
   h_wave->GetXaxis()->SetLabelSize(0.05);
   h_wave->GetXaxis()->SetTitleSize(0.05);
   h_wave->GetYaxis()->SetLabelSize(0.05);
   h_wave->GetYaxis()->SetTitleSize(0.05);
   h_wave->GetYaxis()->SetTitleOffset(1.);
   h_wave->GetYaxis()->SetTitle("ADC counts");
   h_wave->GetXaxis()->SetTitle("time (counts)");
   h_wave->GetXaxis()->SetNdivisions(505);
   h_wave->SetLineColor(2);
    h_wave->Draw("hist");
    
    Double_t hmax = h_wave->GetMaximum();
    
    sprintf (string,"xpos=%.1f cm, R=%.2f\n",xpos,Rup);
    printf("string=%s",string);
    TLatex *t1 = new TLatex(0.2,0.92,string);
    t1->SetNDC();
    t1->SetTextSize(0.05);
    // t1->Draw();
    
    c1->cd(2);
    
   // TF1 *pulse1 = pulse->DrawCopy("same");
    TH1 *h_wave_cumulative = h_wave->GetCumulative();
    h_wave_cumulative->Scale(1./hmax);
    h_wave_cumulative->GetYaxis()->SetTitle("Integral/Peak");
    h_wave_cumulative->GetXaxis()->SetTitle("time (counts)");
    h_wave_cumulative->Draw();
    
    listfile << "bin\tInteg/Peak" << endl;
    
    Int_t ndim = h_wave_cumulative->GetNbinsX();
    
    for(j=0;j<ndim;j++) {
        Double_t xbin = h_wave_cumulative->GetBinCenter(j+1);
        Double_t content = h_wave_cumulative->GetBinContent(j+1);
        printf ("j=%d, xbin=%f, content=%f\n",j,xbin-20.5,content);
        // listfile << xbin << std::setprecision(2) << content << endl;
        listfile << xbin << std::fixed << std::setprecision(2) << "\t" << content << endl;    }
    
    listfile.close();
    
   sprintf (string,"plot_waveform_integral.pdf");
    c1->SaveAs(string);
    sprintf (string,"plot_waveform_integral.png");
    c1->SaveAs(string);
}