void sipm_data_read_pmts(void)
{
// 
// Input the .csv data file from excel file into arrays for processing.
//

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

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
//gROOT->LoadMacro("$ROOTSYS/test/libEvent.so"); 
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 line[1024];
   char filename[132];
   char name[10];
   Double_t pi=3.14159;
   Int_t j,jj;


  // input data for radiometer readings

   # define nps1 29;   
   Double_t f1[nps1]={ 1000000.0,  2000000.0,  3000000.0,  4000000.0,  5000000.0,  6000000.0,  7000000.0,  8000000.0,  9000000.0,  10000000.0,  1000000.0,  500000.0,  800000.0,  700000.0,  900000.0,  600000.0,  400000.0,  300000.0,  200000.0,  100000.0,  90000.0,  80000.0,  70000.0,  60000.0,  50000.0,  40000.0,  30000.0,  20000.0,  10000.0};
   // Double_t radiometer1[nps1];
   Double_t radiopulse1[nps1]={0.000035,  0.000036,  0.000037,  0.000039,  0.000040,  0.000041,  0.000042,  0.000044,  0.000046,  0.000048,  0.000035,  0.000033,  0.000034,  0.000033,  0.000034,  0.000033,  0.000032,  0.000032,  0.000031,  0.000031,  0.000030,  0.000030,  0.000030,  0.000030,  0.000030,  0.000030,  0.000029,  0.000030,  0.000032};

   # define nps2 28;   
   Double_t f2[nps2]={10000.0,  20000.0,  30000.0,  40000.0,  50000.0,  60000.0,  70000.0,  80000.0,  90000.0,  100000.0,  200000.0,  300000.0,  400000.0,  500000.0,  600000.0,  700000.0,  800000.0,  900000.0,  1000000.0,  2000000.0,  3000000.0,  4000000.0,  5000000.0,  6000000.0,  7000000.0,  8000000.0,  9000000.0,  10000000.0};
   Double_t radiopulse2[nps2]={0.000077,  0.000095,  0.000090,  0.000087,  0.000088,  0.000087,  0.000087,  0.000086,  0.000087,  0.000090,  0.000090,  0.000090,  0.000090,  0.000090,  0.000090,  0.000090,  0.000091,  0.000091,  0.000094,  0.000096,  0.000099,  0.000102,  0.000102,  0.000105,  0.000107,  0.000110,  0.000114,  0.000121};

   // normalize radiometer readings to nppe

  // Double_t pulse_norm1=1.71e6;
  // Double_t pulse_norm2=1.71e6/3;
  Double_t pulse_norm1=1e7;
  Double_t pulse_norm2=1e7;
  for (j=0;j<nps1;j++) {
         radiopulse1[j]=radiopulse1[j]*pulse_norm1;
        }

  for (j=0;j<nps2;j++) {
         radiopulse2[j]=radiopulse2[j]*pulse_norm2;
        }

   
   sprintf(filename,"%s","foreltonregularpmts.csv");
   printf ("Read from filename=%s\n",filename);
   ifstream in1 (filename);

   // #define npts 22;     // 1954
   // #define npts 23;     // 1928
   #define npts 19;     // pmts
   // #define npts 20;     // 1977
   // #define npts 20;     // 1998
   Int_t maxlines=1000;
   Double_t gain2[npts];
   Double_t gain3[npts];
   Double_t npho23[npts];
   Double_t tofres1[npts];
   Double_t tofres2[npts];
   Double_t nppe1[npts];
   Double_t nppe2[npts];
   Double_t tr1[npts];
   Double_t tr2[npts];
   Double_t co2[npts];
   Double_t co3[npts];
   Double_t Run[npts];
   Double_t Hz[npts];
   Double_t tofres_ratio[npts];
   Double_t tr_ratio[npts];
   Double_t tofres_ratio2[npts];
   Double_t Inv_co2[npts];
   Double_t Inv_co3[npts];

   Double_t xmin=0;
   Double_t xmax=1.e7;   // zero to 10MHz
   Double_t ymin=0.;
   Double_t ymax=10;


   // read over headers

   in1.getline(line,1024,'\r');
   in1.getline(line,1024,'\r');
   in1.getline(line,1024,'\r');
   in1.getline(line,1024,'\r');

   jj = 0;
          //     while(!in1.eof()) {

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> gain2[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> gain3[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> npho23[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> tofres1[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> tofres2[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> nppe1[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> nppe2[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> tr1[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> tr2[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> co2[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> co3[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> Run[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

          in1.getline(line,1024,'\r');
          jj++;
         printf ("line=%s\n",line);
         stringstream input(line);
         input >> name;
         printf ("first name=%s\n",name);
         for (j=0;j<npts;j++) {
              input >> Hz[j];
              printf ("name=%s, gain2=%f\n",name,gain2[j]);
             }

        // compute derived quantities

       for (j=0;j<npts;j++) {
                     tofres1[j] * = 0.035;  // 35 ps per count
                     tofres2[j] * = 0.035;  // 35 ps per count
                     tofres_ratio[j] = tofres2[j]/tofres1[j];
                     tr_ratio[j] = tr2[j]/tr1[j];
                     tofres_ratio2[j] = tofres_ratio[j]*tofres_ratio[j];
                     Inv_co2[j] = 1/co2[j];
                     Inv_co3[j] = 1/co3[j];
                     printf ("tofres_ratio=%f, tr_ratio=%f, tofres_ratio2=%f, inv_co2=%f, inv_co3=%f\n\n",tofres_ratio[j],tr_ratio[j],tofres_ratio2[j],Inv_co2[j],Inv_co3[j]);
                     }


   //
   TCanvas *c1 = new TCanvas("c1","c1 sipm_data_read_pmts Canvas",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,2);

   c1->cd(1);  
   c1_1->SetGridx();
   c1_1->SetGridy();
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);

   TGraph *gr_tofres_ratio = new TGraph (npts,Hz,tofres_ratio);
   /*t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();*/


  TF1 *pol = new TF1("uncert",pol_func,0,10,4);
  pol->SetLineColor(2);
  Double_t a = 2;
  Double_t b = 0;
  Double_t c = 0;
  Double_t d = 0;
  pol->SetParameters(a,b,c,d);
  pol->SetParNames("a","b","c","d");
  gr_tofres_ratio->Fit(pol);
  

   xmin = 1.e3;
   xmax = 1.e7;
   c1_1->SetLogx();
   ymin = 0;
   ymax = 5;
   gr_tofres_ratio->SetTitle("Rpmts: Bias=73.6V, Amp=2.08V");
   gr_tofres_ratio->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_tofres_ratio->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_tofres_ratio->GetXaxis()->SetTitleSize(0.05);
   gr_tofres_ratio->GetYaxis()->SetTitleSize(0.05);
   gr_tofres_ratio->GetYaxis()->SetTitleOffset(1.5);
   gr_tofres_ratio->GetYaxis()->SetTitle("Resolution Ratio");
   gr_tofres_ratio->GetXaxis()->SetTitle("Frequency (Hz)");
   gr_tofres_ratio->GetXaxis()->SetNdivisions(5);
   gr_tofres_ratio->SetMarkerColor(1);
   gr_tofres_ratio->SetMarkerStyle(21);
   gr_tofres_ratio->SetMarkerSize(0.7);
   gr_tofres_ratio->Draw("AP");

   c1->cd(2); 
   c1_2->SetGridx();
   c1_2->SetGridy();
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   c1_2->SetLogx();

   TGraph *gr_tr_ratio = new TGraph (npts,Hz,tr_ratio);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 5;
   gr_tr_ratio->SetTitle("");
   gr_tr_ratio->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_tr_ratio->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_tr_ratio->GetXaxis()->SetTitleSize(0.05);
   gr_tr_ratio->GetYaxis()->SetTitleSize(0.05);
   gr_tr_ratio->GetYaxis()->SetTitleOffset(1.5);
   gr_tr_ratio->GetYaxis()->SetTitle("Tr Ratio");
   gr_tr_ratio->GetXaxis()->SetTitle("Frequency (Hz)");
   gr_tr_ratio->GetXaxis()->SetNdivisions(5);
   gr_tr_ratio->SetMarkerColor(1);
   gr_tr_ratio->SetMarkerStyle(21);
   gr_tr_ratio->SetMarkerSize(0.7);
   gr_tr_ratio->Draw("AP");


   c1->cd(3); 
   c1_3->SetGridx();
   c1_3->SetGridy();
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   c1_3->SetLogx();

   TGraph *gr_inv_co2 = new TGraph (npts,Hz,Inv_co2);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 0.2;
   gr_inv_co2->SetTitle("");
   gr_inv_co2->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_inv_co2->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_inv_co2->GetXaxis()->SetTitleSize(0.05);
   gr_inv_co2->GetYaxis()->SetTitleSize(0.05);
   gr_inv_co2->GetYaxis()->SetTitleOffset(1.5);
   gr_inv_co2->GetYaxis()->SetTitle("1/co2");
   gr_inv_co2->GetXaxis()->SetTitle("Frequency (Hz)");
   gr_inv_co2->GetXaxis()->SetNdivisions(5);
   gr_inv_co2->SetMarkerColor(1);
   gr_inv_co2->SetMarkerStyle(21);
   gr_inv_co2->SetMarkerSize(0.7);
   gr_inv_co2->Draw("AP");


   c1->cd(4); 
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);
   c1_4->SetLogx();

   TGraph *gr_inv_co3 = new TGraph (npts,Hz,Inv_co3);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 0.2;
   gr_inv_co3->SetTitle("");
   gr_inv_co3->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_inv_co3->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_inv_co3->GetXaxis()->SetTitleSize(0.05);
   gr_inv_co3->GetYaxis()->SetTitleSize(0.05);
   gr_inv_co3->GetYaxis()->SetTitleOffset(1.5);
   gr_inv_co3->GetYaxis()->SetTitle("1/co3");
   gr_inv_co3->GetXaxis()->SetTitle("Frequency (Hz)");
   gr_inv_co3->GetXaxis()->SetNdivisions(5);
   gr_inv_co3->SetMarkerColor(1);
   gr_inv_co3->SetMarkerStyle(21);
   gr_inv_co3->SetMarkerSize(0.7);
   gr_inv_co3->Draw("AP");


   //
   TCanvas *c2 = new TCanvas("c2","c2 sipm_data_read_pmts Canvas",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,2);

   c2->cd(1);  
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   c2_1->SetLogx();

   TGraph *gr_npho23 = new TGraph (npts,Hz,npho23);
   /*t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();*/

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 2000;
   gr_npho23->SetTitle("Rpmts: Bias=73.6V, Amp=2.08V");
   gr_npho23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_npho23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_npho23->GetXaxis()->SetTitleSize(0.05);
   gr_npho23->GetYaxis()->SetTitleSize(0.05);
   gr_npho23->GetYaxis()->SetTitleOffset(1.5);
   gr_npho23->GetYaxis()->SetTitle("Npho 23");
   gr_npho23->GetXaxis()->SetTitle("Frequency (Hz)");
   gr_npho23->GetXaxis()->SetNdivisions(5);
   gr_npho23->SetMarkerColor(1);
   gr_npho23->SetMarkerStyle(21);
   gr_npho23->SetMarkerSize(0.7);
   gr_npho23->Draw("AP");

   c2->cd(2); 
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   // c2_2->SetLogx();

   TGraph *gr_ratio_n23 = new TGraph (npts,npho23,tofres_ratio);

   xmin = 0;
   xmax = 2000;
   ymin = 0;
   ymax = 5;
   gr_ratio_n23->SetTitle("");
   gr_ratio_n23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_ratio_n23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_ratio_n23->GetXaxis()->SetTitleSize(0.05);
   gr_ratio_n23->GetYaxis()->SetTitleSize(0.05);
   gr_ratio_n23->GetYaxis()->SetTitleOffset(1.5);
   gr_ratio_n23->GetYaxis()->SetTitle("Resolution Ratio");
   gr_ratio_n23->GetXaxis()->SetTitle("Npho23");
   gr_ratio_n23->GetXaxis()->SetNdivisions(5);
   gr_ratio_n23->SetMarkerColor(1);
   gr_ratio_n23->SetMarkerStyle(21);
   gr_ratio_n23->SetMarkerSize(0.7);
   gr_ratio_n23->Draw("AP");


   c2->cd(3); 
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   c2_3->SetLogy();
   c2_3->SetLogx();

   TGraph *gr_res1_vs_co2 = new TGraph (npts,Inv_co2,tofres1);
   TGraph *gr_res2_vs_co3 = new TGraph (npts,Inv_co3,tofres2);

   TLegend *leg1 = new TLegend(0.5,0.8,0.85,0.9);
   leg1->AddEntry(gr_res1_vs_co2,"tofres1","p");
   leg1->AddEntry(gr_res2_vs_co3,"tofres2","p");


   xmin = 0.04;
   xmax = 0.2;
   ymin = 0.05;
   ymax = 0.2;
   Double_t dummyx[3]={0.01,0.1,0.3};
   Double_t dummyy[3]={0,0.1,0.2};
   TGraph *dum = new TGraph(3,dummyx,dummyy);
   dum->SetTitle("");
   dum->GetXaxis()->SetRangeUser(xmin,xmax);
   dum->GetYaxis()->SetRangeUser(ymin,ymax);
   dum->GetXaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleOffset(1.5);
   dum->GetYaxis()->SetTitle("Resolution (TDC counts)");
   dum->GetXaxis()->SetTitle("1/co2, 1/co3");
   dum->Draw("AP");
   

   xmin = 0.01;
   xmax = 0.2;
   ymin = 0.04;
   ymax = 0.4;
   gr_res2_vs_co3->SetTitle("");
   gr_res2_vs_co3->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_res2_vs_co3->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_res2_vs_co3->GetXaxis()->SetTitleSize(0.05);
   gr_res2_vs_co3->GetYaxis()->SetTitleSize(0.05);
   gr_res2_vs_co3->GetYaxis()->SetTitleOffset(1.5);
   gr_res2_vs_co3->GetYaxis()->SetTitle("Resolution (TDC counts)");
   gr_res2_vs_co3->GetXaxis()->SetTitle("1/co2, 1/co3");
   // gr_res2_vs_co3->GetXaxis()->SetNdivisions(5);
   gr_res2_vs_co3->SetMarkerColor(4);
   gr_res2_vs_co3->SetMarkerStyle(21);
   gr_res2_vs_co3->SetMarkerSize(0.7);
   gr_res2_vs_co3->Draw("SameP");
   gr_res1_vs_co2->SetMarkerColor(2);
   gr_res1_vs_co2->SetMarkerStyle(21);
   gr_res1_vs_co2->SetMarkerSize(0.7);
   gr_res1_vs_co2->Draw("SameP");
   leg1->Draw();


   c2->cd(4); 
   c2_4->SetGridx();
   c2_4->SetGridy();
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   c2_4->SetLogy();
   c2_4->SetLogx();

   TGraph *gr_res1_vs_hz = new TGraph (npts,Hz,tofres1);
   TGraph *gr_res2_vs_hz = new TGraph (npts,Hz,tofres2);

   TLegend *leg2 = new TLegend(0.5,0.8,0.85,0.9);
   leg2->AddEntry(gr_res1_vs_hz,"R2083","p");
   leg2->AddEntry(gr_res2_vs_hz,"R9779","p");


   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0.04;
   ymax = 0.4;
   gr_res1_vs_hz->SetTitle("");
   gr_res1_vs_hz->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_res1_vs_hz->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_res1_vs_hz->GetXaxis()->SetTitleSize(0.05);
   gr_res1_vs_hz->GetYaxis()->SetTitleSize(0.05);
   gr_res1_vs_hz->GetYaxis()->SetTitleOffset(1.5);
   gr_res1_vs_hz->GetYaxis()->SetTitle("Resolution (ns)");
   gr_res1_vs_hz->GetXaxis()->SetTitle("Rate (Hz)");
   // gr_res1_vs_hz->GetXaxis()->SetNdivisions(5);
   gr_res1_vs_hz->GetYaxis()->SetMoreLogLabels();
   gr_res1_vs_hz->SetMarkerColor(2);
   gr_res1_vs_hz->SetMarkerStyle(21);
   gr_res1_vs_hz->SetMarkerSize(0.7);
   gr_res1_vs_hz->Draw("AP");
   leg2->Draw();

   gr_res2_vs_hz->SetTitle("");
   gr_res2_vs_hz->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_res2_vs_hz->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_res2_vs_hz->GetXaxis()->SetTitleSize(0.05);
   gr_res2_vs_hz->GetYaxis()->SetTitleSize(0.05);
   gr_res2_vs_hz->GetYaxis()->SetTitleOffset(1.5);
   gr_res2_vs_hz->GetYaxis()->SetTitle("Resolution (ns)");
   gr_res2_vs_hz->GetXaxis()->SetTitle("Rate (Hz)");
   // gr_res2_vs_hz->GetXaxis()->SetNdivisions(5);
   gr_res2_vs_hz->GetYaxis()->SetMoreLogLabels();
   gr_res2_vs_hz->SetMarkerColor(4);
   gr_res2_vs_hz->SetMarkerStyle(21);
   gr_res2_vs_hz->SetMarkerSize(0.7);
   gr_res2_vs_hz->Draw("SameP");
   leg2->Draw();


   //
   TCanvas *c3 = new TCanvas("c3","c3 sipm_data_read_pmts Canvas",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   c3->Divide(2,2);

   c3->cd(1);  
   c3_1->SetGridx();
   c3_1->SetGridy();
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);
   c3_1->SetLogx();
   // c3_1->SetLogy();

   TGraph *gr_res1_vs_npho23 = new TGraph (npts,npho23,tofres1);
   TGraph *gr_res2_vs_npho23 = new TGraph (npts,npho23,tofres2);

   TLegend *leg3 = new TLegend(0.15,0.8,0.45,0.9);
   leg3->AddEntry(gr_res1_vs_hz,"tofres1","p");
   leg3->AddEntry(gr_res2_vs_hz,"tofres2","p");

   xmin = 100;
   xmax = 2000;
   ymin = 1;
   ymax = 10;
   gr_res1_vs_npho23->SetTitle("Rpmts: Bias=73.6V, Amp=2.08V");
   gr_res1_vs_npho23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_res1_vs_npho23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_res1_vs_npho23->GetXaxis()->SetTitleSize(0.05);
   gr_res1_vs_npho23->GetYaxis()->SetTitleSize(0.05);
   gr_res1_vs_npho23->GetYaxis()->SetTitleOffset(1.5);
   gr_res1_vs_npho23->GetYaxis()->SetTitle("Resolution (TDC counts)");
   gr_res1_vs_npho23->GetXaxis()->SetTitle("Npho23");
   // gr_res1_vs_npho23->GetYaxis()->SetNdivisions(9);
   // gr_res1_vs_npho23->GetXaxis()->SetNdivisions(10);
   gr_res1_vs_npho23->SetMarkerColor(2);
   gr_res1_vs_npho23->SetMarkerStyle(21);
   gr_res1_vs_npho23->SetMarkerSize(0.7);
   gr_res1_vs_npho23->GetXaxis()->SetMoreLogLabels();
   gr_res1_vs_npho23->Draw("AP");
   leg3->Draw();

  TF1 *sqrt = new TF1("uncert",sqrt_func,0,2000,2);
  sqrt->SetLineColor(2);
  Double_t a = 80;
  Double_t b = 0.6;
  sqrt->SetParameters(a,b);
  sqrt->FixParameter(1,b);
  sqrt->SetParNames("a","b");
  fitmin=500;
  fitmax=2000;
  gr_res1_vs_npho23->Fit(sqrt,"","",fitmin,fitmax);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 0.3;

   Double_t dummyx[3]={1.e3,1.e4,1.e7};
   Double_t dummyy[3]={-1,-1,-1};
   TGraph *dum = new TGraph(3,dummyx,dummyy);
   dum->SetTitle("Pmts: R9779 -HV1500, R2083 -HV1900");
   dum->GetXaxis()->SetRangeUser(xmin,xmax);
   dum->GetYaxis()->SetRangeUser(ymin,ymax);
   dum->GetXaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleOffset(1.5);
   dum->GetYaxis()->SetTitle("Resolution (ns)");
   dum->GetXaxis()->SetTitle("Rate (Hz)");
   // dum->GetXaxis()->SetMoreLogLabels();
   dum->Draw("AP");

  gr_res1_vs_hz->Draw("SameP");
  gr_res2_vs_hz->Draw("SameP");
  leg2->Draw();

   c3->cd(2);  
   c3_2->SetGridx();
   c3_2->SetGridy();
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   c3_2->SetLogx();
   // c3_2->SetLogy();

   gr_res2_vs_npho23->SetTitle("");
   gr_res2_vs_npho23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_res2_vs_npho23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_res2_vs_npho23->GetXaxis()->SetTitleSize(0.05);
   gr_res2_vs_npho23->GetYaxis()->SetTitleSize(0.05);
   gr_res2_vs_npho23->GetYaxis()->SetTitleOffset(1.5);
   gr_res2_vs_npho23->GetYaxis()->SetTitle("Resolution (TDC counts)");
   gr_res2_vs_npho23->GetXaxis()->SetTitle("Npho23");
   // gr_res2_vs_npho23->GetXaxis()->SetNdivisions(5);
   gr_res2_vs_npho23->SetMarkerColor(4);
   gr_res2_vs_npho23->SetMarkerStyle(21);
   gr_res2_vs_npho23->SetMarkerSize(0.7);
   gr_res2_vs_npho23->GetXaxis()->SetMoreLogLabels();
   // gr_res2_vs_npho23->Draw("AP");
   // leg3->Draw();


  TF1 *sqrt = new TF1("uncert",sqrt_func,0,2000,2);
  sqrt->SetLineColor(4);
  Double_t a = 80;
  Double_t b = 0.6;
  sqrt->SetParameters(a,b);
  sqrt->FixParameter(1,b);
  sqrt->SetParNames("a","b");
  fitmin=500;
  fitmax=2000;
  gr_res2_vs_npho23->Fit(sqrt,"","",fitmin,fitmax);


   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   // ymax = 0.3;
   ymax = 2000;

   Double_t dummyx[3]={1.e3,1.e4,1.e7};
   Double_t dummyy[3]={-1,-1,-1};
   TGraph *dum = new TGraph(3,dummyx,dummyy);
   dum->SetTitle("");
   dum->GetXaxis()->SetRangeUser(xmin,xmax);
   dum->GetYaxis()->SetRangeUser(ymin,ymax);
   dum->GetXaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleOffset(1.5);
   dum->GetYaxis()->SetTitle("Number of photoelectrons");
   dum->GetXaxis()->SetTitle("Rate (Hz)");
   // dum->GetXaxis()->SetMoreLogLabels();
   dum->Draw("AP");

  gr_npho23->Draw("Samep");
  leg2->Draw();


   c3->cd(3); 
   c3_3->SetGridx();
   c3_3->SetGridy();
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);
   c3_3->SetLogx();

   TGraph *gr_ratio_hz = new TGraph (npts,Hz,tofres_ratio);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 3;

   Double_t dummyx[3]={1.e3,1.e4,1.e7};
   Double_t dummyy[3]={-1,-1,-1};
   TGraph *dum = new TGraph(3,dummyx,dummyy);
   dum->SetTitle("");
   dum->GetXaxis()->SetRangeUser(xmin,xmax);
   dum->GetYaxis()->SetRangeUser(ymin,ymax);
   dum->GetXaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleSize(0.05);
   dum->GetYaxis()->SetTitleOffset(1.5);
   dum->GetYaxis()->SetTitle("Resolution Ratio");
   dum->GetXaxis()->SetTitle("Rate (Hz)");
   dum->Draw("AP");

   gr_ratio_hz->SetMarkerColor(1);
   gr_ratio_hz->SetMarkerStyle(21);
   gr_ratio_hz->SetMarkerSize(0.7);
   gr_ratio_hz->Draw("SameP");

  TF1 *sratio = new TF1("sratio",sratio_func,0,1.e7,3);
  sratio->SetLineColor(4);
  Double_t Npixel = 1000;
  Double_t tp = 10.e-6;
  Double_t rise_ratio=1;
  sratio->SetParameters(Npixel,tp,rise_ratio);
  // sratio->FixParameter(1,b);
  sratio->SetParNames("Nphotons","tp","Rtau");
  fitmin=1000;
  fitmax=1e7;
  gr_ratio_hz->Fit(sratio,"","",fitmin,fitmax);

  c3->cd(4);
   c3_4->SetGridx();
   c3_4->SetGridy();
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   c3_4->SetLogx();

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 2000;

   TGraph *pulse1 = new TGraph (nps1,f1,radiopulse1);
   pulse1->SetTitle("Radiometer");
   pulse1->GetXaxis()->SetRangeUser(xmin,xmax);
   pulse1->GetYaxis()->SetRangeUser(ymin,ymax);
   pulse1->GetXaxis()->SetTitleSize(0.05);
   pulse1->GetYaxis()->SetTitleSize(0.05);
   pulse1->GetYaxis()->SetTitleOffset(1.5);
   pulse1->GetYaxis()->SetTitle("fW/cm2/pulse");
   pulse1->GetXaxis()->SetTitle("Frequency (Hz)");
   pulse1->GetXaxis()->SetNdivisions(5);
   pulse1->SetMarkerColor(2);
   pulse1->SetMarkerStyle(20);
   pulse1->SetMarkerSize(0.7);
   pulse1->Draw("Ap");

   TGraph *pulse2 = new TGraph (nps2,f2,radiopulse2);
   pulse2->SetMarkerColor(2);
   pulse2->SetMarkerStyle(24);
   pulse2->SetMarkerSize(0.7);
   pulse2->Draw("Samep");


  c1->SaveAs("sipm_data_read_pmtsb_c1.pdf");
  c2->SaveAs("sipm_data_read_pmtsb_c2.pdf");
  c3->SaveAs("sipm_data_read_pmtsb_c3.pdf");


  // close input file
}
Double_t pol_func (Double_t *x, Double_t *par) {

     // return uncertainty in y0 = x0*a + b,  value as computed using covariance matrix elements

     Double_t a=par[0];
     Double_t b=par[1];
     Double_t c=par[2];
     Double_t d=par[2];
     Double_t x1=x[0]; 

     char string[256];

     Double_t func=0;
    
     func = a + b*x1 + c*x1*x1 + d*x1*x1*x1;

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

     // return  y = a/sqrt(x) + b, with b=0.6 TDC channels (1 TDC channel= 35ps)

     Double_t a=par[0];
     Double_t b=par[1];
     Double_t x1=x[0]; 

     char string[256];

     Double_t func=0;
    
     if (x1 > 0) {
          func = a /sqrt(x1) + b;
          }

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

     // return uncertainty in y0 = x0*a + b,  value as computed using covariance matrix elements

     Double_t N=par[0];    // number of photons hitting the PMT
     Double_t tp=par[1];   // dead time parameter of SiPM
     Double_t rise_ratio=par[2];   // ratio of rise time parameters ts/tp =  SiPM/reference
     Double_t Qp=0.17;    // QE of reference pmt at 470 nm
     Double_t Qs=0.21;    // QE of SiPM at 470 nm
     Double_t fsipm=1.0;    // geometric area of sipm compact (0.48 compact, 0.59 production) relative to the PMT
     Double_t Ntot=3600*16;   // Total number of pixes in the SiPM array
     Double_t p0 = 566.2;    // parameter of Npe vs f
     Double_t p1 = 0.0001857;  // parameter of Npe vs f
     // Double_t p1 = 0.1857;  // parameter of Npe vs f
     Double_t f=x[0];

     // Function assumed is Eq. 13 of GlueX-doc-1807
     // Parameterization of the ratio of the SiPM resolution to the reference timing resolution
    

     char string[256];

     Double_t func=0;
     // Int_t exponent = f*tp;
     Float_t exponent = f*tp;

     Double_t Nsipm = fsipm*N;
     // Nsipm = fsipm*N*(1+p1*f/p0);
     // Double_t R = 1;   // remove dead time dependent term
     Double_t R = pow((1- (Qs*Nsipm/Ntot)),exponent);
     func = sqrt(Ntot/(Qp*N)) * sqrt(1-exp(-Qs*Nsipm/Ntot))*sqrt(R)/rise_ratio;

     // printf("f=%g exponent=%g Nsipm=%g tp=%g Rtau=%g R=%g func=%g\n",f,exponent,Nsipm,tp,rise_ratio,R,func);

     return 1./func;
}