void plot_wrapping (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 npts1 5;
   #define npts2 4;
   #define npts3 6;
   #define jpts 10


   Double_t ynone[npts1]={-1,427,553,450,-1};
   Double_t xnone[npts1]={-1,1,1,1,5};
   Double_t ytedlar[npts2]={-1,491,492,5};
   Double_t xtedlar[npts2]={-1,1.95,2.05,-1};
   Double_t ypaper[npts2]={-1,528,504,5};
   Double_t xpaper[npts2]={-1,3,3,-1};

   Double_t ynone2[npts1]={-1,210,192,224,-1};
   Double_t xnone2[npts1]={-1,1,1,1,5};
   Double_t ytedlar2[npts1]={-1,180,167,215,-1};
   Double_t xtedlar2[npts1]={-1,2,2,2,-1};
   Double_t ypaper2[npts1]={-1,171,177,224,-1};
   Double_t xpaper2[npts1]={-1,3,3,3,-1};

   Double_t ynone3[npts2]={-1,54.5,46.5,-1};
   Double_t xnone3[npts2]={-1,1,1,-1};
   Double_t ytedlar3[npts2]={-1,44.6,49.3,-1};
   Double_t xtedlar3[npts2]={-1,2,2,-1};
   Double_t ypaper3[npts2]={-1,48.2,51.6,-1};
   Double_t xpaper3[npts2]={-1,3,3,-1};

   Double_t ynone4[npts1]={-1,44.5,47.9,48.5,-1};
   Double_t xnone4[npts1]={-1,1,1,1,6};
   Double_t ytedlar4[npts1]={-1,41.7,43.3,45.2,-1};
   Double_t xtedlar4[npts1]={-1,2,2,2,6};
   Double_t ypaper4[npts1]={-1,48.1,47.1,41.3,-1};
   Double_t xpaper4[npts1]={-1,3,3,3,6};
   Double_t ymylar4[npts1]={-1,44,49.2,45.3,-1};
   Double_t xmylar4[npts1]={-1,4,4,4,6};

   Double_t factor1=2.5;
   Double_t factor2=10;

   for (j=0;j<npts1;j++) {
          ynone2[j] *= factor1;
          ytedlar2[j] *= factor1;
          ypaper2[j] *= factor1;
          }

   for (j=0;j<npts2;j++) {
          ynone3[j] *= factor2;
          ytedlar3[j] *= factor2;
          ypaper3[j] *= factor2;
          }

   // compute mean and sigma

  Double_t S1none=0;
  Double_t S2none=0;
  Double_t S1tedlar=0;
  Double_t S2tedlar=0;
  Double_t S1paper=0;
  Double_t S2paper=0;
  Double_t S1mylar=0;
  Double_t S2mylar=0;
  Double_t sum=0;

  for (j=1;j<npts1-1;j++) {
         S1none += ynone4[j];
         S2none += ynone4[j]*ynone4[j];
         S1tedlar += ytedlar4[j];
         S2tedlar += ytedlar4[j]*ytedlar4[j];
         S1paper += ypaper4[j];
         S2paper += ypaper4[j]*ypaper4[j];
         S1mylar += ymylar4[j];
         S2mylar += ymylar4[j]*ymylar4[j];
         sum += 1;
         }
    Double_t mean_none = S1none/sum;
    Double_t sigma_none= sqrt((S2none - sum*mean_none*mean_none)/(sum-1));
    Double_t mean_tedlar = S1tedlar/sum;
    Double_t sigma_tedlar= sqrt((S2tedlar - sum*mean_tedlar*mean_tedlar)/(sum-1));
    Double_t mean_paper = S1paper/sum;
    Double_t sigma_paper= sqrt((S2paper - sum*mean_paper*mean_paper)/(sum-1));
    Double_t mean_mylar = S1mylar/sum;
    Double_t sigma_mylar= sqrt((S2mylar - sum*mean_mylar*mean_mylar)/(sum-1));
    cout << " mean none " << mean_none << " sigma_none " << sigma_none << endl;
    cout << " mean tedlar " << mean_tedlar << " sigma_tedlar " << sigma_tedlar << endl;
    cout << " mean paper " << mean_paper << " sigma_paper " << sigma_paper << endl;
    cout << " mean mylar " << mean_mylar << " sigma_mylar " << sigma_mylar << endl;

   Double_t means[4];
   Double_t sigmas[4];
   means[0] = mean_none;
   means[1] = mean_tedlar;
   means[2] = mean_paper;
   means[3] = mean_mylar;
   sigmas[0] = sigma_none;
   sigmas[1] = sigma_tedlar;
   sigmas[2] = sigma_paper;
   sigmas[3] = sigma_mylar; 
   Double_t xmeans[4]={1.1,2.1,3.1,4.1};
   Double_t xsigmas[4]={0.01,0.01,0.01,0.01};

   TGraph *mean = new TGraphErrors(4,xmeans,means,xsigmas,sigmas);


   TGraph *none= new TGraph(npts1,xnone,ynone); 
   TGraph *tedlar= new TGraph(npts2,xtedlar,ytedlar);   
   TGraph *paper= new TGraph(npts2,xpaper,ypaper);   

   TGraph *none2= new TGraph(npts1,xnone2,ynone2); 
   TGraph *tedlar2= new TGraph(npts1,xtedlar2,ytedlar2);   
   TGraph *paper2= new TGraph(npts1,xpaper2,ypaper2);   

   TGraph *none3= new TGraph(npts2,xnone3,ynone3); 
   TGraph *tedlar3= new TGraph(npts2,xtedlar3,ytedlar3);   
   TGraph *paper3= new TGraph(npts2,xpaper3,ypaper3);     

   TGraph *none4= new TGraph(npts1,xnone4,ynone4); 
   TGraph *tedlar4= new TGraph(npts1,xtedlar4,ytedlar4);   
   TGraph *paper4= new TGraph(npts1,xpaper4,ypaper4);     
   TGraph *mylar4= new TGraph(npts1,xmylar4,ymylar4);   

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

   TLegend *leg = new TLegend(0.15,0.70,0.50,0.95);

   Double_t xmin=0;
   Double_t xmax=4;
   Double_t ymin=0;
   Double_t ymax=800;
   none->SetTitle("");
   none->GetXaxis()->SetRangeUser(xmin,xmax);
   none->GetYaxis()->SetRangeUser(ymin,ymax);
   none->GetXaxis()->SetTitleSize(0.04);
   none->GetYaxis()->SetTitleSize(0.04);
   none->GetYaxis()->SetTitleOffset(1.5);
   // none->GetYaxis()->SetTitleOffset(2.0);
   none->GetXaxis()->SetTitle("Wrapping Type");
   // none->GetYaxis()->SetTitle("Amplitude (mV)");
   none->GetYaxis()->SetTitle("Normalized Amplitude (mV)");
   none->GetXaxis()->SetNdivisions(505);
   none->SetMarkerColor(1);
   none->SetMarkerStyle(21);

   leg->AddEntry(none,"No Wrapping (22 mm)","p");
   none->SetMarkerColor(4);
   none->SetMarkerStyle(21);
   none->Draw("Ap");

   leg->AddEntry(tedlar,"Tedlar  (22 mm)","p");
   tedlar->SetMarkerColor(4);
   tedlar->SetMarkerStyle(21);
   tedlar->Draw("Samep");

   leg->AddEntry(paper,"Acid-free Paper (22 mm)","p");
   paper->SetMarkerColor(4);
   paper->SetMarkerStyle(21);
   paper->Draw("Samep");

   leg->AddEntry(none2,"No Wrapping  x 2.5 (80 mm)","p");
   none2->SetMarkerColor(1);
   none2->SetMarkerStyle(25);
   none2->Draw("Samep");

   leg->AddEntry(tedlar2,"Tedlar  x 2.5 (80 mm)","p");
   tedlar2->SetMarkerColor(1);
   tedlar2->SetMarkerStyle(25);
   tedlar2->Draw("Samep");

   leg->AddEntry(paper2,"Acid-free Paper  x 2.5 (80 mm)","p");
   paper2->SetMarkerColor(1);
   paper2->SetMarkerStyle(25);
   paper2->Draw("Samep");

   leg->AddEntry(none3,"No Wrapping x10 (diffuser)","p");
   none3->SetMarkerColor(2);
   none3->SetMarkerStyle(20);
   none3->Draw("Samep");

   leg->AddEntry(tedlar3,"Tedlar x10 (diffuser)","p");
   tedlar3->SetMarkerColor(2);
   tedlar3->SetMarkerStyle(20);
   tedlar3->Draw("Samep");

   leg->AddEntry(paper3,"Acid-free Paper x10 (diffuser)","p");
   paper3->SetMarkerColor(2);
   paper3->SetMarkerStyle(20);
   paper3->Draw("Samep");

   leg->Draw();

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

   TLegend *leg1 = new TLegend(0.15,0.75,0.50,0.90);

   Double_t xmin=0;
   Double_t xmax=5;
   Double_t ymin=0;
   Double_t ymax=80;
   none4->SetTitle("");
   none4->GetXaxis()->SetRangeUser(xmin,xmax);
   none4->GetYaxis()->SetRangeUser(ymin,ymax);
   none4->GetXaxis()->SetTitleSize(0.04);
   none4->GetYaxis()->SetTitleSize(0.04);
   none4->GetYaxis()->SetTitleOffset(1.5);
   // none4->GetYaxis()->SetTitleOffset(2.0);
   none4->GetXaxis()->SetTitle("Wrapping Type");
   none4->GetYaxis()->SetTitle("Amplitude (mV)");
   none4->GetXaxis()->SetNdivisions(510);
   none4->SetMarkerColor(1);
   none4->SetMarkerStyle(21);

   leg1->AddEntry(none4,"No Wrapping (diffuser)","p");
   none4->SetMarkerColor(1);
   none4->SetMarkerStyle(21);
   none4->Draw("Ap");

   leg1->AddEntry(tedlar4,"Tedlar  (diffuser)","p");
   tedlar4->SetMarkerColor(2);
   tedlar4->SetMarkerStyle(21);
   tedlar4->Draw("Samep");

   leg1->AddEntry(paper4,"Acid-free Paper (diffuser)","p");
   paper4->SetMarkerColor(4);
   paper4->SetMarkerStyle(21);
   paper4->Draw("Samep");

   leg1->AddEntry(mylar4," Al Mylar (diffuser)","p");
   mylar4->SetMarkerColor(3);
   mylar4->SetMarkerStyle(21);
   mylar4->Draw("Samep");

   mean->SetMarkerColor(2);
   mean->SetMarkerStyle(20);
   mean->Draw("Samep");

   Double_t barline= (mean_none+mean_tedlar+mean_paper+mean_mylar)/4;
   TLine *bar = new TLine(xmin,barline,xmax,barline);
    bar->SetLineColor(1);
    bar->Draw();

   leg1->Draw();

   sprintf(filename,"plot_wrapping_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"plot_wrapping_c1.png");
   c1->SaveAs(filename);
   sprintf(filename,"plot_wrapping_c2.eps");
   c2->SaveAs(filename);
   sprintf(filename,"plot_wrapping_c2.png");
   c2->SaveAs(filename);
   /*sprintf(filename,"plot_wrapping_c3.eps");
   c3->SaveAs(filename);
   sprintf(filename,"plot_wrapping_c3.png");
   c3->SaveAs(filename);*/

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

   // Fit ln(thermal conductivity) as a funciton of ln(Temperature)

   Double_t p0=par[0];
   Double_t p1=par[1];
   Double_t p2=par[2];
   Double_t p3=par[3];
   Double_t p4=par[4];

   Double_t x1=log(x[0]);
   Double_t pi=3.14159;

   char string[256];

   Double_t func;

   func = p0 + p1*x1 + p2*x1*x1 + p3*x1*x1*x1 + p4*x1*x1*x1*x1;

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

   if (func < 50) {
       return exp(func);
      }
   else {
      return exp(50);
     }
}