// #include <TMath.h>
#include <TRandom.h>
#include <MakeTSpline.C>
void fit_bcal_pulse(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//


gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetOptStat(111111);
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;

	// These points are taken by mapping out figure 6
// of GlueX-doc-1795 (the "5ns" rise time pulse).

#define npts 32

Double_t time[npts] ={
-7.00,
-4.92,
-0.75,
4.24,
9.37,	// 5
9.71,
10.05,
10.39,
10.73,
11.07,
12.03,
12.92,
14.08,
16.05,	// 10
18.56,
21.22,
24.31,
28.44,
32.92,	// 15
38.64,
42.07,
47.10,
52.59,
58.22,	// 20
67.95,
71.74,
79.85,
120.0,
125.0,	// 25
130.0,
160.0,
180.0
};


Double_t mV[npts] ={
0,
0,
0,
0,
-1,	// 5
-1.2,
-1.4,
-1.6,
-1.8,
-2,
-49.2,
-159.0,
-375.9,
-711.4,	// 10
-902.2,
-989.0,
-1020.8,
-960.1,
-850.2,	// 15
-694.0,
-601.5,
-539.3,
-419.9,
-300.5,	// 20
-142.5,
-100.6,
-61.6,
-8.0,
-4.0,	// 25
-2.0,
-1.0,
0.0
};


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

   Double_t xmin=-20;
   Double_t xmax=180;
   // use fitted values
   Double_t pulse_mu=0.058;  // result of fit  - good recovery to baseline
   Double_t pulse_sigma=5.91;     // result of fit  - slower rise time
   // Double_t pulse_mu=0.04;  //   eyeball fit - baseline recovers too late.
   // Double_t pulse_sigma=3.5;     // eyeball fit - matches fast rise time  
   Double_t gain=-34970.;   

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

   Double_t sum=0;
   Double_t delta=4;

   //
   TCanvas *c1 = new TCanvas("c1","c1 fit_bcal_pulse",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);*/
   pulse->SetTitle("");
   // pulse->GetYaxis()->SetRangeUser(0,2*pedestal);   
   pulse->GetXaxis()->SetTitleSize(0.04);
   pulse->GetYaxis()->SetTitleSize(0.04);
   pulse->GetYaxis()->SetTitleOffset(1.5);
   pulse->GetXaxis()->SetTitle("Time (ns)");
   pulse->GetYaxis()->SetTitle("Arb. Units");
   pulse->GetXaxis()->SetNdivisions(505);
   pulse->Draw();


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

   TGraph *laser_data = new TGraph(npts,time,mV);
   laser_data->SetMarkerStyle(20);
   laser_data->GetXaxis()->SetNdivisions(505);
   laser_data->GetXaxis()->SetTitle("Time (ns)");
   laser_data->Draw("Ap");
   pulse->Draw("Csame");
   // laser_data->Fit(pulse);

   // get spline fit

  TSpline *fspline = MakeTSpline();
  fspline->Draw("Csame");


   sprintf(filename,"fit_bcal_pulse_c2.pdf");
   c2->SaveAs(filename);
   sprintf(filename,"fit_bcal_pulse_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]-18;

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