void fit_pde(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];
   Int_t j,jj;
   
   #define npts 9;

   // input data set first point (actually -1) to -10 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 V[npts]={0,1,1.5,2,2.5,3,3.5,4,5};
   Double_t xerr[npts]={npts*0.1};

   // Dark rates vs temperature at fixed overvoltage.

   Double_t pde[npts]={-1,6,11.5,17,19,20.5,23,25,-1};
   Double_t dr_err[npts]={npts*0.5};

   Double_t xmin=0;
   Double_t xmax=5;
   Double_t ymin=0;
   Double_t ymax=30;

   //
   TCanvas *c1 = new TCanvas("c1","PDE vs Over-voltage",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   // c1->SetLogy();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);


   TGraph *plot = new TGraph (npts,V,pde);

   TLegend *leg = new TLegend(0.15,0.85,0.75,0.9);

   plot->SetTitle("");
   plot->GetXaxis()->SetRangeUser(xmin,xmax);
   plot->GetYaxis()->SetRangeUser(ymin,ymax);
   plot->GetXaxis()->SetTitleSize(0.04);
   plot->GetYaxis()->SetTitleSize(0.04);
   plot->GetYaxis()->SetTitleOffset(1.5);
   plot->GetXaxis()->SetTitle("Over-voltage (V)");
   plot->GetYaxis()->SetTitle("PDE (%)");
   plot->GetXaxis()->SetNdivisions(6);

   plot->SetMarkerColor(2);
   plot->SetMarkerStyle(20);
   plot->SetMarkerSize(1.5);
   plot->Draw("AP");

   Double_t xminfit = 1;
   Double_t xmaxfit = 3.5;
   plot->Fit("pol3","","",xminfit,xmaxfit);
   Double_t a = pol3->GetParameter(0);
   Double_t b = pol3->GetParameter(1);
   Double_t c = pol3->GetParameter(2);
   Double_t d = pol3->GetParameter(3);
   sprintf (string,"Lightfoot 460 nm SensL 1mm\n");
   leg->AddEntry(plot,string,"p");
   leg->Draw();

   sprintf (string,"a=%.2f\n",a);
   printf ("string=%s\n",string);
   t1 = new TLatex(0.2,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"b=%.2f\n",b);
   printf ("string=%s\n",string);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"c=%.2f\n",c);
   printf ("string=%s\n",string);
   t1 = new TLatex(0.2,0.7,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"d=%.2f\n",d);
   printf ("string=%s\n",string);
   t1 = new TLatex(0.2,0.65,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   //
   TCanvas *c2 = new TCanvas("c2","PDE vs Over-voltage",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   // c2->SetLogy();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   Double_t slope[npts];
   slope[0] = -1;
   slope[npts-1] = -1;
   for (j=1;j<npts-1;j++) {

      // slope % change per mV

      slope[j] = pol3->Derivative(V[j])*100/(pde[j]*1000);
      printf ("j=%d V=%f, slope=%f\n",j,V[j],slope[j]);
      }

   TGraph *sl = new TGraph(npts,V,slope);
   ymin = 0;
   ymax = 0.1;

   sl->SetTitle("");
   sl->GetXaxis()->SetRangeUser(xmin,xmax);
   sl->GetYaxis()->SetRangeUser(ymin,ymax);
   sl->GetXaxis()->SetTitleSize(0.04);
   sl->GetYaxis()->SetTitleSize(0.04);
   sl->GetYaxis()->SetTitleOffset(2);
   sl->GetXaxis()->SetTitle("Over-voltage (V)");
   sl->GetYaxis()->SetTitle("PDE Change (%/mV)");
   sl->GetXaxis()->SetNdivisions(6);

   sl->SetMarkerColor(2);
   sl->SetMarkerStyle(20);
   sl->SetMarkerSize(1.5);
   sl->Draw("AP");
   leg->Draw();

   //

   c1->SaveAs("fit_pde.eps");
   c1->SaveAs("fit_pde.gif");
   c2->SaveAs("fit_pde_slope.eps");
   c2->SaveAs("fit_pde_slope.gif");

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

   Double_t a=par[0];
   Double_t b=par[1];
   // Double_t c=par[2];
   Double_t c=0;
   Double_t T0=20;
   Double_t x1=x[0]-T0;

   char string[256];

   Double_t exp1 = a*exp(b*x1 + c*x1*x1);
   
   /*sprintf (string,"exp1=%f a=%f b=%f c=%f\n",exp1,a,b,c);
   printf ("string=%s",string);*/

   return exp1;
}