void plot_MeVee (void)
{
// 
// plot the results of light transmission for various light guide wrappings
//

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

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];
   char filename[80];
   Int_t j,jj;
   
   #define npts 12;
   Double_t xmin=0;
   Double_t xmax=100;
   Double_t ymin=0;
   Double_t ymax=100;

   Double_t proton[npts]={-1,1,2,3,10,15,16.6,20,40,50,100,120};
   Double_t pout[npts]={-1,0.2,0.55,0.9,5.2,8.7,10,12.5,28,34,70,-1};

   TGraph *output = new TGraph(npts,proton,pout);
   TF1 *mevee = new TF1("mevee",mevee_func,xmin,xmax,1);
    mevee->SetParameter(0,1);

   TF1 *parab = new TF1("parab",parab_func,xmin,xmax,2);
    parab->SetParameter(0.7,0);

   // output->Fit("pol2","","",0,100);
   /* output->Fit("mevee","","",0,20);
   Double_t alpha = mevee->GetParameter(0);
   cout << "alpha=" << alpha << endl;*/

   output->Fit("parab","","",0,100);
   Double_t alpha = parab->GetParameter(0);
   Double_t p1 = parab->GetParameter(1);
   Double_t alpha_err = parab->GetParError(0);
   Double_t p1_err = parab->GetParError(1);
    cout << "alpha=" << alpha << " alpha_err=" << alpha_err << " p1=" << p1 << " p1_err=" << p1_err << endl;

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

   TLegend *leg = new TLegend(0.15,0.85,0.60,0.90);

   output->SetTitle("");
   output->GetXaxis()->SetRangeUser(xmin,xmax);
   output->GetYaxis()->SetRangeUser(ymin,ymax);
   output->GetXaxis()->SetTitleSize(0.04);
   output->GetYaxis()->SetTitleSize(0.04);
   output->GetYaxis()->SetTitleOffset(1.5);
   // output->GetYaxis()->SetTitleOffset(2.0);
   output->GetXaxis()->SetTitle("Proton Energy (MeV)");
   output->GetYaxis()->SetTitle("Light Output (MeV_{ee})");
   output->GetXaxis()->SetNdivisions(505);
   output->SetMarkerColor(1);
   output->SetMarkerStyle(20);

   leg->AddEntry(output,"G.V O'Rielly NIM A368 (1996) 745","p");
   output->SetMarkerColor(4);
   output->SetMarkerStyle(20);
   output->Draw("Ap");
   // mevee->Draw("same");
   parab->Draw("same");

   leg->Draw();

   sprintf (string,"f(x) = (%.2f #pm %.2f)x + (%.4f #pm %.4f)x^{2}\n",alpha,alpha_err,p1,p1_err);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.69,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();

   sprintf(filename,"plot_MeVee_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"plot_MeVee_c1.png");
   c1->SaveAs(filename);

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

   // Parameterization of the light output for protons.

   Double_t alpha=par[0];
   // Double_t p1=par[1];
   // Double_t p2=par[2];

   Double_t x1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   func = 0.5*(exp(alpha*x1)+exp(-alpha*x1));
   func = log(func);

   /*sprintf (string,"x1=%f func=%f\n",x1,func);
   printf ("string=%s",string);*/

   return func;
}
Double_t parab_func (Double_t *x, Double_t *par) {

   // Parameterization of the light output for protons.

   Double_t alpha=par[0];
   Double_t p1=par[1];
   // Double_t p2=par[2];

   Double_t x1=x[0];
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   func = alpha*x1 + p1*x1*x1;

   /*sprintf (string,"x1=%f func=%f\n",x1,func);
   printf ("string=%s",string);*/

   return func;
}