void sipm_data_read_302(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.71e7;
  Double_t pulse_norm2=1.71e7/3;
  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","9-20-2011_302.csv");
   printf ("Read from filename=%s\n",filename);
   ifstream in1 (filename);

   # define npts 23;        // 302
   // #define npts 22;     // 1954
   // #define npts 23;     // 1928
   // #define npts 19;     // 1928
   // #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 sigma1[npts];   // TDC units are 35 ps/count
   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 tpnpho23[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

       jj = 0;
       // Double_t Hzcut=1.6e6;
       Double_t Hzcut=10.e6;
       for (j=0;j<npts;j++) {
                     // reduce hits based on rate
                     if (Hz[j] < Hzcut) {
                        gain2[jj] = gain2[j];
                        gain3[jj] = gain3[j];
                        npho23[jj] = npho23[j];
                        tofres1[jj] = tofres1[j];
                        sigma1[jj] = tofres1[j]*0.035;   // convert TDC counts to ns.
                        tofres2[jj] = tofres2[j];
                        nppe1[jj] = nppe1[j];
                        nppe2[jj] = nppe2[j];
                        tr1[jj] = tr1[j];
                        tr2[jj] = tr2[j];
                        co2[jj] = co2[j];
                        co3[jj] = co3[j];
                        Run[jj] = Run[j];
                        Hz[jj] = Hz[j];

                        tofres_ratio[jj] = tofres2[jj]/tofres1[jj];
                         tr_ratio[jj] = tr2[jj]/tr1[jj];
                        tofres_ratio2[jj] = tofres_ratio[jj]*tofres_ratio[jj];
                        Inv_co2[jj] = 1/co2[jj];
                        Inv_co3[jj] = 1/co3[jj];
                        tpnpho23[jj]=tofres1[jj]*sqrt(npho23[jj]);
                        printf ("tofres_ratio=%f, tr_ratio=%f, tofres_ratio2=%f, inv_co2=%f, inv_co3=%f\n\n",tofres_ratio[jj],tr_ratio[jj],tofres_ratio2[jj],Inv_co2[jj],Inv_co3[jj]);
                        jj++;
                       }
                     }
    Int_t npts1 = jj;
    printf ("\nSelect data npts1=%d\n",npts1);


   //
   TCanvas *c1 = new TCanvas("c1","c1 sipm_data_read_302 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 (npts1,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("R302: Bias=71.6V, -HV2800");
   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 (npts1,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 (npts1,Hz,Inv_co2);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 1;
   // gr_inv_co2->SetTitle("");
   gr_inv_co2->SetTitle("R302: Bias=71.6V, -HV2800");
   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("Amplitude (arb. units)");
   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 (npts1,Hz,Inv_co3);

   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 0;
   ymax = 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_302 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 (npts1,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("R302: Bias=71.6V, -HV2800");
   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("Npe");
   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 (npts1,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 (npts1,Inv_co2,tofres1);
   TGraph *gr_res2_vs_co3 = new TGraph (npts1,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.2;
   xmax = 2;
   ymin = 1;
   ymax = 40;
   Double_t dummyx[3]={0.01,0.1,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 = 1;
   ymax = 40;
   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 (npts1,Hz,tofres1);
   TGraph *gr_res2_vs_hz = new TGraph (npts1,Hz,tofres2);

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


   xmin = 1.e3;
   xmax = 1.e7;
   ymin = 1;
   ymax = 40;
   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 (TDC counts)");
   gr_res1_vs_hz->GetXaxis()->SetTitle("Hz");
   // gr_res1_vs_hz->GetXaxis()->SetNdivisions(5);
   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");
   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_302 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 (npts1,npho23,tofres1);
   TGraph *gr_res2_vs_npho23 = new TGraph (npts1,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("R302: Bias=71.6V, -HV2800");
   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("Npe");
   // 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);

   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);


   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 (npts1,Hz,tofres_ratio);

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


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

   c4->cd(1);  
   c4_1->SetGridx();
   c4_1->SetGridy();
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);
   c4_1->SetLogx();
   // c4_1->SetLogy();

   gr_inv_co2->SetMarkerColor(4);
   gr_inv_co2->Draw("AP");


   c4->cd(2);  
   c4_2->SetGridx();
   c4_2->SetGridy();
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   c4_2->SetLogx();
   // c4_2->SetLogy();



   gr_npho23->Draw("AP");
   gr_npho23->SetMarkerColor(4);
  fitmin=1000;
  fitmax=6.e6;
  gr_npho23->Fit("pol1","","",fitmin,fitmax);

   TGraph *pulse1 = new TGraph (nps1,f1,radiopulse1);
   pulse1->SetMarkerColor(2);
   pulse1->SetMarkerStyle(20);
   pulse1->SetMarkerSize(0.7);
   pulse1->Draw("Samep");

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

   c4->cd(3); 
   c4_3->SetGridx();
   c4_3->SetGridy();
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);
   c4_3->SetLogy();
   c4_3->SetLogx();

   TGraph *gr_sigma1_vs_npho23 = new TGraph (npts1,npho23,sigma1);

   xmin = 100;
   xmax = 2000;
   ymin = 0.05;
   ymax = 0.2;
   gr_sigma1_vs_npho23->SetTitle("R302: Bias=71.6V, -HV2800");
   gr_sigma1_vs_npho23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_sigma1_vs_npho23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_sigma1_vs_npho23->GetXaxis()->SetTitleSize(0.05);
   gr_sigma1_vs_npho23->GetYaxis()->SetTitleSize(0.05);
   gr_sigma1_vs_npho23->GetYaxis()->SetTitleOffset(1.5);
   gr_sigma1_vs_npho23->GetYaxis()->SetTitle("Resolution (ns)");
   gr_sigma1_vs_npho23->GetXaxis()->SetTitle("Npe");
   // gr_sigma1_vs_npho23->GetYaxis()->SetNdivisions(9);
   // gr_sigma1_vs_npho23->GetXaxis()->SetNdivisions(10);
   gr_sigma1_vs_npho23->SetMarkerColor(2);
   gr_sigma1_vs_npho23->SetMarkerStyle(21);
   gr_sigma1_vs_npho23->SetMarkerSize(0.7);
   gr_sigma1_vs_npho23->GetXaxis()->SetMoreLogLabels();
   gr_sigma1_vs_npho23->Draw("AP");

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

   c4->cd(4); 
   c4_4->SetGridx();
   c4_4->SetGridy();
   c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);
   c4_4->SetLogy();
   c4_4->SetLogx();


   TGraph *gr_tpnpho23_vs_npho23 = new TGraph (npts1,npho23,tpnpho23);

   xmin = 100;
   xmax = 2000;
   ymin = 20;
   ymax = 200;
   gr_tpnpho23_vs_npho23->SetTitle("R302: Bias=71.6V, -HV2800");
   gr_tpnpho23_vs_npho23->GetXaxis()->SetRangeUser(xmin,xmax);
   gr_tpnpho23_vs_npho23->GetYaxis()->SetRangeUser(ymin,ymax);
   gr_tpnpho23_vs_npho23->GetXaxis()->SetTitleSize(0.05);
   gr_tpnpho23_vs_npho23->GetYaxis()->SetTitleSize(0.05);
   gr_tpnpho23_vs_npho23->GetYaxis()->SetTitleOffset(1.5);
   gr_tpnpho23_vs_npho23->GetYaxis()->SetTitle("Resolution x Sqrt(Npe)");
   gr_tpnpho23_vs_npho23->GetXaxis()->SetTitle("Npe");
   // gr_tpnpho23_vs_npho23->GetYaxis()->SetNdivisions(9);
   // gr_tpnpho23_vs_npho23->GetXaxis()->SetNdivisions(10);
   gr_tpnpho23_vs_npho23->SetMarkerColor(2);
   gr_tpnpho23_vs_npho23->SetMarkerStyle(21);
   gr_tpnpho23_vs_npho23->SetMarkerSize(0.7);
   gr_tpnpho23_vs_npho23->GetXaxis()->SetMoreLogLabels();
   gr_tpnpho23_vs_npho23->Draw("AP");


  c1->SaveAs("sipm_data_read_302a_c1.pdf");
  c2->SaveAs("sipm_data_read_302a_c2.pdf");
  c3->SaveAs("sipm_data_read_302a_c3.pdf");
  c4->SaveAs("sipm_data_read_302a_c4.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;
}