void atten_waveform(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;

   // 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;

  Double_t xmin=0;
  Double_t xmax=16;
  Double_t ymin=0;
  Double_t ymax=0.2;
  
  TH1F *hslope = new TH1F("hslope","Waveform slope", nbins,ymin,ymax);
  TH1F *hsigma = new TH1F("hsigma","Waveform sigma", nbins,ymin,ymax);

   Int_t npts=16;
   Double_t mu_data[16]={0.079,0.083,0.073,0.079,0.100,0.076,0.076,0.084,0.094,0.067,0.073,0.089,0.078,0.077,0.071,0.081};
   Double_t mu_err[16]={0.008,0.008,0.007,0.003,0.008,0.004,0.009,0.005,0.009,0.006,0.004,0.006,0.009,0.010,0.007,0.007};
   Double_t sigma_data[16]={1.2,0.85,0.82,0.89,0.88,0.79,0.94,0.77,1.17,0.05,0.76,0.94,1.01,1.85,1.12,1.05};
   Double_t sigma_err[16]={0.30,0.21,0.27,0.15,0.23,0.19,0.38,0.15,0.42,0.60,0.12,0.21,0.36,0.59,0.34,0.33};
   Double_t ndx_data[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
   Double_t ndx_err[16]={npts*0.001};

   // compute weighted means


   Double_t wm=0;
   Double_t w=0;

   for (jj=0;jj<npts;jj++) {
          Double_t wi = 1/(mu_err[jj]*mu_err[jj]);
          wm += mu_data[jj]*wi;
          w += wi;
         }
    Double_t mu_mean=wm/w;

   wm = 0;
   w = 0;
   for (jj=0;jj<npts;jj++) {
          Double_t wi = 1/(sigma_err[jj]*sigma_err[jj]);
          wm += sigma_data[jj]*wi;
          w += wi;
         }
    Double_t mu_sigma=wm/w;

   printf(" mu_mean=%f, mu_sigma=%f\n",mu_mean,mu_sigma);
          

   TGraphErrors *gmu = new TGraphErrors(npts,ndx_data,mu_data,ndx_err,mu_err);
   TGraphErrors *gsigma = new TGraphErrors(npts,ndx_data,sigma_data,ndx_err,sigma_err);


   TCanvas *c0 = new TCanvas("c0", "c0",200,10,700,700);
   c0->SetBorderMode(0);
   c0->SetFillColor(0);

   c0->Divide(1,2);
   c0->cd(1);
   c0_1->SetBorderMode(0);
   c0_1->SetFillColor(0);


   gmu->GetXaxis()->SetRangeUser(xmin,xmax);
   gmu->GetYaxis()->SetRangeUser(ymin,ymax);
   gmu->GetXaxis()->SetTitleSize(0.07);
   gmu->GetYaxis()->SetTitleSize(0.07);
   gmu->GetXaxis()->SetTitle("Channel");
   gmu->GetYaxis()->SetTitle("Slope (1/counts)");
   gmu->SetMarkerColor(2);
   gmu->SetMarkerStyle(20);
   gmu->SetTitle("");
   gmu->Draw("AP");
   gmu->Fit("pol0");

   TLine *line = new TLine(xmin,mu_mean,xmax,mu_mean);
    line->SetLineColor(2);
    line->SetLineWidth(2);
    line->Draw();
   sprintf (string,"mean=%.4f\n",mu_mean);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   c0->cd(2);
   c0_2->SetBorderMode(0);
   c0_2->SetFillColor(0);

  Double_t ymin=0;
  Double_t ymax=2;

   gsigma->GetXaxis()->SetRangeUser(xmin,xmax);
   gsigma->GetYaxis()->SetRangeUser(ymin,ymax);
   gsigma->GetXaxis()->SetTitleSize(0.07);
   gsigma->GetYaxis()->SetTitleSize(0.07);
   gsigma->GetXaxis()->SetTitle("Channel");
   gsigma->GetYaxis()->SetTitle("sigma (counts)");
   gsigma->SetMarkerColor(2);
   gsigma->SetMarkerStyle(20);
   gsigma->SetTitle("");
   gsigma->Draw("AP");
   gsigma->Fit("pol0");

   TLine *line = new TLine(xmin,mu_sigma,xmax,mu_sigma);
    line->SetLineColor(2);
    line->SetLineWidth(2);
    line->Draw();
   sprintf (string,"mean=%.3f\n",mu_sigma);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


   // define histogram to hold pulses

  Double_t xmin=-100;
  Double_t xmax=300;

  TH1F *hup_wave= new TH1F("hup_wave","Up waveform", nbins,xmin,xmax);
  TH1F *hupR_wave= new TH1F("hupR_wave","Up Reflected waveform", nbins,xmin,xmax);
  TH1F  *hdn_wave= new TH1F("hdn_wave","Down waveform", nbins,xmin,xmax);
  TH1F  *hdnR_wave= new TH1F("hdnR_wave","Down Reflected waveform", nbins,xmin,xmax);
  TH1F *hupSum_wave= new TH1F("hupSum_wave","Up Sum waveform", nbins,xmin,xmax);
  TH1F *hdnSum_wave= new TH1F("hdnSum_wave","Down Sum waveform", nbins,xmin,xmax);

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

   Double_t xpmin=0;
   Double_t xpmax=100;

   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 up_time = -xpos/veff + time0;
         Double_t up_atten = exp(-(L/2+xpos)/lambda);
         Double_t up_x = up_time/4 + t0;
         Double_t up_y = up_atten*(pulse->Eval(up_x) - pedestal);
         hup_wave->Fill(time0,up_y);

         Double_t upR_time = -(-xpos+L)/veff + time0 ;
         Double_t upR_atten = exp(-(L/2-xpos+L)/lambda);
         Double_t upR_x = upR_time/4 + t0;
         Double_t upR_y = upR_atten*Rup*(pulse->Eval(upR_x) - pedestal);
         hupR_wave->Fill(time0,upR_y);

         Double_t dn_time = xpos/veff + time0;
         Double_t dn_atten = exp(-(L/2-xpos)/lambda);
         Double_t dn_x = dn_time/4 + t0;
         Double_t dn_y = dn_atten*(pulse->Eval(dn_x) - pedestal);
         hdn_wave->Fill(time0,dn_y);

         Double_t dnR_time = -(xpos+L)/veff + time0;
         Double_t dnR_atten = exp(-(L/2+xpos+L)/lambda);
         Double_t dnR_x = dnR_time/4 + t0;
         Double_t dnR_y = dnR_atten*Rdn*(pulse->Eval(dnR_x) - pedestal);
         hdnR_wave->Fill(time0,dnR_y);

         hupSum_wave->Fill(time0,up_y);
         hupSum_wave->Fill(time0,upR_y);
         hdnSum_wave->Fill(time0,dn_y);
         hdnSum_wave->Fill(time0,dnR_y);

         // printf (" jj=%d time0=%f up_time=%f up_x=%f up_y=%f dn_time=%f dn_x=%f dn_y=%f\n",jj,time0,up_time,up_x,up_y,dn_time,dn_x,dn_y);
         // printf (" jj=%d time0=%f upR_time=%f upR_x=%f upR_y=%f dnR_time=%f dnR_x=%f dnR_y=%f\n",jj,time0,upR_time,upR_x,upR_y,dnR_time,dnR_x,dnR_y);
         }
         
  
   TCanvas *c1 = new TCanvas("c1", "c1",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();

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

   hup_wave->SetTitle("");
   hup_wave->GetXaxis()->SetLabelSize(0.05);
   hup_wave->GetXaxis()->SetTitleSize(0.05);
   hup_wave->GetYaxis()->SetLabelSize(0.05);
   hup_wave->GetYaxis()->SetTitleSize(0.05);
   hup_wave->GetYaxis()->SetTitleOffset(1.5);
   hup_wave->GetYaxis()->SetTitle("UP ADC counts");
   hup_wave->GetXaxis()->SetTitle("time (ns)");
   hup_wave->GetXaxis()->SetNdivisions(505);
   hup_wave->SetLineColor(2);
   hup_wave->Draw();   

   sprintf (string,"xpos=%.1f cm, R=%.2f\n",xpos,Rup);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);

   hdn_wave->SetTitle("");
   hdn_wave->GetXaxis()->SetLabelSize(0.05);
   hdn_wave->GetXaxis()->SetTitleSize(0.05);
   hdn_wave->GetYaxis()->SetLabelSize(0.05);
   hdn_wave->GetYaxis()->SetTitleSize(0.05);
   hdn_wave->GetYaxis()->SetTitleOffset(1.5);
   hdn_wave->GetYaxis()->SetTitle("DOWN ADC counts");
   hdn_wave->GetXaxis()->SetTitle("time (ns)");
   hdn_wave->GetXaxis()->SetNdivisions(505);
   hdn_wave->SetLineColor(4);
   hdn_wave->Draw();

   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);

   hupR_wave->SetTitle("");
   hupR_wave->GetXaxis()->SetLabelSize(0.05);
   hupR_wave->GetXaxis()->SetTitleSize(0.05);
   hupR_wave->GetYaxis()->SetLabelSize(0.05);
   hupR_wave->GetYaxis()->SetTitleSize(0.05);
   hupR_wave->GetYaxis()->SetTitleOffset(1.5);
   hupR_wave->GetYaxis()->SetTitle("UP ADC counts");
   hupR_wave->GetXaxis()->SetTitle("time (ns)");
   hupR_wave->GetXaxis()->SetNdivisions(505);
   hupR_wave->SetLineColor(2);
   hupR_wave->Draw();

   c1->cd(4);
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);

   hdnR_wave->SetTitle("");
   hdnR_wave->GetXaxis()->SetLabelSize(0.05);
   hdnR_wave->GetXaxis()->SetTitleSize(0.05);
   hdnR_wave->GetYaxis()->SetLabelSize(0.05);
   hdnR_wave->GetYaxis()->SetTitleSize(0.05);
   hdnR_wave->GetYaxis()->SetTitleOffset(1.5);
   hdnR_wave->GetYaxis()->SetTitle("DOWN ADC counts");
   hdnR_wave->GetXaxis()->SetTitle("time (ns)");
   hdnR_wave->GetXaxis()->SetNdivisions(505);
   hdnR_wave->SetLineColor(4);
   hdnR_wave->Draw();   

   TCanvas *c2 = new TCanvas("c2", "c2",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   c2->SetGridy();

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

   hupSum_wave->SetTitle("");
   hupSum_wave->GetXaxis()->SetLabelSize(0.05);
   hupSum_wave->GetXaxis()->SetTitleSize(0.05);
   hupSum_wave->GetYaxis()->SetLabelSize(0.05);
   hupSum_wave->GetYaxis()->SetTitleSize(0.05);
   hupSum_wave->GetYaxis()->SetTitleOffset(1.5);
   hupSum_wave->GetYaxis()->SetTitle("UP ADC counts");
   hupSum_wave->GetXaxis()->SetTitle("time (ns)");
   hupSum_wave->GetXaxis()->SetNdivisions(505);
   hupSum_wave->SetLineColor(2);
   hupSum_wave->Draw();   

   sprintf (string,"xpos=%.1f cm, R=%.2f\n",xpos,Rup);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);

   hdnSum_wave->SetTitle("");
   hdnSum_wave->GetXaxis()->SetLabelSize(0.05);
   hdnSum_wave->GetXaxis()->SetTitleSize(0.05);
   hdnSum_wave->GetYaxis()->SetLabelSize(0.05);
   hdnSum_wave->GetYaxis()->SetTitleSize(0.05);
   hdnSum_wave->GetYaxis()->SetTitleOffset(1.5);
   hdnSum_wave->GetYaxis()->SetTitle("DOWN ADC counts");
   hdnSum_wave->GetXaxis()->SetTitle("time (ns)");
   hdnSum_wave->GetXaxis()->SetNdivisions(505);
   hdnSum_wave->SetLineColor(4);
   hdnSum_wave->Draw();


   sprintf (string,"atten_waveform_%.0f.pdf(",xpos);
   c0->SaveAs(string);
   sprintf (string,"atten_waveform_%.0f.pdf",xpos);
   c1->SaveAs(string);
   sprintf (string,"atten_waveform_%.0f.pdf)",xpos);
   c2->SaveAs(string);

}
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 if {

         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;
}