void plot_pde (void)
{
// 
// plot the pde from measurements in Eckard NIM A620 (2010) 217, Fig 8.
//

// #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 mpts 19;
   #define npts 20;
   #define jpts 10


   Double_t Vover[mpts]={-1,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2};
   Double_t pde_v[mpts]={-1,9,11.5,14,15.5,17.5,19.2,21,22.5,24,25,26.5,27.5,28.5,29.5,30.5,31.,32.5,-1};

   Double_t lambda[npts]={300,360,380,400,420,440,460,480,500,560,600,660,700,760,800,860,900,960,1000,1500};
   Double_t pde_w[npts]={-1,28,28.4,30.1,31,32.5,32,31,30.1,24.3,20.9,15.8,12.8,9.0,7.5,4.8,3.0,1.6,0.9,-1};

   // plot pde vs overvoltage

   TCanvas *c1 = new TCanvas("c1","c1 plot_pde",200,10,700,700);
 
   TGraph *pde_vover = new TGraph (mpts,Vover,pde_v);

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

   Double_t xmin=-1;
   Double_t xmax=2.5;
   Double_t ymin=0;
   Double_t ymax=50;
   pde_vover->SetTitle("");
   pde_vover->GetXaxis()->SetRangeUser(xmin,xmax);
   pde_vover->GetYaxis()->SetRangeUser(ymin,ymax);
   pde_vover->GetXaxis()->SetTitleSize(0.04);
   pde_vover->GetYaxis()->SetTitleSize(0.04);
   pde_vover->GetYaxis()->SetTitleOffset(1.5);
   // pde_vover->GetYaxis()->SetTitleOffset(2.0);
   pde_vover->GetXaxis()->SetTitle("Over Voltage (V)");
   pde_vover->GetYaxis()->SetTitle("PDE (%)");
   pde_vover->GetXaxis()->SetNdivisions(505);
   pde_vover->SetMarkerColor(1);
   pde_vover->SetMarkerStyle(21);
   pde_vover->Draw("Ap");

   Double_t fitmin=1.0;
   Double_t fitmax=2.;
   pde_vover->Fit("pol1","","",fitmin,fitmax);
   TF1 *fitv1 = pde_vover->GetFunction("pol1");
   Float_t p0 = fitv1->GetParameter(0); 
   Float_t p1 = fitv1->GetParameter(1);
   Float_t p0err = fitv1->GetParError(0); 
   Float_t p1err = fitv1->GetParError(1); 
   Float_t p2,p2err;

   TF1 *fitv1full = new TF1("fitvfull","pol1",-1,fitmax);
   fitvfull->SetParameters(p0,p1);
   fitvfull->Draw("same");
   fitvfull->SetLineColor(4);
   fitvfull->SetLineStyle(2);

   printf ("Fit1 parameters: p0 %f +/-%f p1 %f  +/-%f p2 %f +/-%f\n",p0,p0err,p1,p1err,p2,p2err);

   Double_t fitmin=0;
   Double_t fitmax=2.1;
   pde_vover->Fit("pol2","+","",fitmin,fitmax);
   TF1 *fitv = pde_vover->GetFunction("pol2");
   Float_t p0 = fitv->GetParameter(0); 
   Float_t p1 = fitv->GetParameter(1); 
   Float_t p2 = fitv->GetParameter(2); 
   Float_t p0err = fitv->GetParError(0); 
   Float_t p1err = fitv->GetParError(1); 
   Float_t p2err = fitv->GetParError(2); 
   fitv->SetLineColor(2);

   printf ("Fit parameters: p0 %f +/-%f p1 %f  +/-%f p2 %f +/-%f\n",p0,p0err,p1,p1err,p2,p2err);   

   sprintf (string,"Fig. 8 from Eckert NIMA620 (2010) 217\n");
   printf("string=%s",string);
   t1 = new TLatex(0.1,0.95,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   Double_t R1 = fitv->Eval(0.9)/fitv->Eval(1.4);
   printf ("Ratio 0.9/1.4nm %f\n",R1);
   Double_t R2 = fitv->Eval(0.9)/fitv->Eval(1.1);
   printf ("Ratio 0.9/1.1 nm %f\n",R2);
   Double_t R3 = fitv->Eval(0.9)/fitv->Eval(1.0);
   printf ("Ratio 0.9/1.1 nm %f\n",R3);

   sprintf (string,"R (0.9/1.4) =%.2f\n",R1);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"R (0.9/1.1)=%.2f\n",R2);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"R (0.9/1.0)=%.2f\n",R3);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();


   // plot pde vs wavelength

   TCanvas *c2 = new TCanvas("c2","c2 plot_pde",200,10,700,700);
 
   TGraph *pde_wave = new TGraph (npts,lambda,pde_w);

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

   Double_t xmin=300;
   Double_t xmax=1000;
   Double_t ymin=0;
   Double_t ymax=50;
   pde_wave->SetTitle("");
   pde_wave->GetXaxis()->SetRangeUser(xmin,xmax);
   pde_wave->GetYaxis()->SetRangeUser(ymin,ymax);
   pde_wave->GetXaxis()->SetTitleSize(0.04);
   pde_wave->GetYaxis()->SetTitleSize(0.04);
   pde_wave->GetYaxis()->SetTitleOffset(1.5);
   // pde_wave->GetYaxis()->SetTitleOffset(2.0);
   pde_wave->GetXaxis()->SetTitle("Wavelength (nm)");
   pde_wave->GetYaxis()->SetTitle("PDE (%)");
   pde_wave->GetXaxis()->SetNdivisions(505);
   pde_wave->SetMarkerColor(2);
   pde_wave->SetMarkerStyle(21);
   pde_wave->Draw("Ap");

   Double_t fitmin=360;
   Double_t fitmax=800;
   TF1 *fitw = new TF1("fitw","[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))+[3]",fitmin,fitmax);
   fitw->SetParameters(31,430,207,0);
   pde_wave->Fit(fitw,"","",fitmin,fitmax);
   // pde_wave->Fit("pol3","","",fitmin,fitmax);
   // TF1 *fitw = pde_wave->GetFunction("pol3");
   Float_t p0 = fitw->GetParameter(0); 
   Float_t p1 = fitw->GetParameter(1); 
   Float_t p2 = fitw->GetParameter(2); 
   Float_t p3 = fitw->GetParameter(3); 
   Float_t p0err = fitw->GetParError(0); 
   Float_t p1err = fitw->GetParError(1); 
   Float_t p2err = fitw->GetParError(2);  
   Float_t p3err = fitw->GetParError(3); 

   printf ("Fit parameters: p0 %f +/-%f p1 %f  +/-%f p2 %f +/-%f p3 %f +/-%f\n",p0,p0err,p1,p1err,p2,p2err,p3,p3err);

   sprintf (string,"Fig. 8 from Eckert NIMA620 (2010) 217\n");
   printf("string=%s",string);
   t1 = new TLatex(0.1,0.95,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   Double_t R0 = fitw->Eval(490)/fitw->Eval(465);
   printf ("Ratio 490/470 nm %f\n",R1);
   Double_t R1 = fitw->Eval(490)/fitw->Eval(470);
   printf ("Ratio 490/470 nm %f\n",R1);
   Double_t R2 = fitw->Eval(490)/fitw->Eval(515);
   printf ("Ratio 490/515 nm %f\n",R2);

   sprintf (string,"R (490/465) =%.2f\n",R0);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"R (490/470) =%.2f\n",R1);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"R (490/515)=%.2f\n",R2);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();


   // plot pde vs wavelength

   TCanvas *c3 = new TCanvas("c3","c3 plot_pde",200,10,700,700);
   Double_t dummyx[3]={0,5,10};
   Double_t dummyy[3]={-1,-1,-1};
   
   TMarker *eckert = new TMarker(1,18*fitw->Eval(490)/fitw->Eval(465),22);
   TMarker *dinu = new TMarker(2,27*fitv->Eval(0.9)/fitv->Eval(1.1),21);
   TMarker *musienko = new TMarker(3,18*fitw->Eval(490)/fitw->Eval(515),29);
   TMarker *musienko2 = new TMarker(3,28*(fitw->Eval(490)/fitw->Eval(515))*(fitv->Eval(0.9)/fitv->Eval(1.0)),30);
   TMarker *bonanno1 = new TMarker(5,20*fitv->Eval(0.9)/fitv->Eval(1.4),24);
   TMarker *bonanno2 = new TMarker(5,45*fitv->Eval(0.9)/fitv->Eval(1.4),24);
   TMarker *bonanno3 = new TMarker(5,27*fitv->Eval(0.9)/fitv->Eval(1.4),20);
   TMarker *usm = new TMarker(6,35,26);
   TMarker *jlab = new TMarker(7,19*fitw->Eval(490)/fitw->Eval(470),25);
   TMarker *hamamatsu = new TMarker(8,38.4,25);

   TGraph *compare = new TGraph(3,dummyx,dummyy);   
   TLegend *leg = new TLegend(0.15,0.60,0.40,0.9);

   
   Double_t xmin=0;
   Double_t xmax=10;
   Double_t ymin=0;
   Double_t ymax=50;
   compare->SetTitle("");
   compare->GetXaxis()->SetRangeUser(xmin,xmax);
   compare->GetYaxis()->SetRangeUser(ymin,ymax);
   compare->GetXaxis()->SetTitleSize(0.04);
   compare->GetYaxis()->SetTitleSize(0.04);
   compare->GetYaxis()->SetTitleOffset(1.5);
   // compare->GetYaxis()->SetTitleOffset(2.0);
   compare->GetXaxis()->SetTitle("Reference");
   compare->GetYaxis()->SetTitle("PDE (%)");
   compare->GetXaxis()->SetNdivisions(505);
   compare->SetMarkerColor(2);
   compare->SetMarkerStyle(20);
   compare->Draw("Ap");

   eckert->SetMarkerColor(2);
   eckert->SetMarkerSize(1.5);
   leg->AddEntry(eckert,"Eckert counting","p");
   eckert->Draw();

   dinu->SetMarkerColor(2);
   dinu->SetMarkerSize(1.5);
   leg->AddEntry(dinu,"Dinu counting","p");
   dinu->Draw();

   musienko->SetMarkerColor(2);
   musienko->SetMarkerSize(1.5);
   leg->AddEntry(musienko,"Musienko counting","p");
   musienko->Draw();

   musienko2->SetMarkerColor(2);
   musienko2->SetMarkerSize(1.5);
   leg->AddEntry(musienko2,"Musienko2 counting","p");
   musienko2->Draw();

   bonanno1->SetMarkerColor(2);
   bonanno1->SetMarkerSize(1.5);
   leg->AddEntry(bonanno1,"Bonanno 100U counting","p");
   bonanno1->Draw();

   bonanno2->SetMarkerColor(4);
   bonanno2->SetMarkerSize(1.5);
   leg->AddEntry(bonanno2,"Bonanno 100U current","p");
   bonanno2->Draw();

   bonanno3->SetMarkerColor(4);
   bonanno3->SetMarkerSize(1.5);
   leg->AddEntry(bonanno3,"Bonanno 50U current","p");
   bonanno3->Draw();

   usm->SetMarkerColor(2);
   usm->SetMarkerSize(1.5);
   leg->AddEntry(usm,"USM #92 counting","p");
   usm->Draw();

   jlab->SetMarkerColor(2);
   jlab->SetMarkerSize(1.5);
   leg->AddEntry(jlab,"JLab #92 counting","p");
   jlab->Draw();

   hamamatsu->SetMarkerColor(4);
   hamamatsu->SetMarkerSize(1.5);
   leg->AddEntry(hamamatsu,"Hamamatsu #92  current","p");
   hamamatsu->Draw();   
 
   leg->Draw();

   sprintf(filename,"plot_pde_c1.pdf");
   c1->SaveAs(filename);
   sprintf(filename,"plot_pde_c1.png");
   c1->SaveAs(filename);
   sprintf(filename,"plot_pde_c2.pdf");
   c2->SaveAs(filename);
   sprintf(filename,"plot_pde_c2.png");
   c2->SaveAs(filename);
   sprintf(filename,"plot_pde_c3.pdf");
   c3->SaveAs(filename);
   sprintf(filename,"plot_pde_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);
     }
}