Double_t Nonlin_int_func (Double_t *x, Double_t *par) {
    
    // Returns the pulse height V=x1 corrected non-linear correction for fininte number of SiPM pixels.
    
    Double_t k = par[0];            // number of pixels/count
    Double_t RAreaPeak = par[1];            // Ratio of area to peak
    Double_t Npixels = par[2];      // number of pixels in one SiPM array
    Double_t layer = par[3];        // layer number for the correction
    Double_t k2 = par[4];           // scale the non-linearity
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    if (x1 <= 0) return 0;
    
    Npixels *= layer;
    Double_t Ntrue = x1;
    
    Double_t Nmeas = Npixels*(1-exp(-k2*Ntrue/Npixels));
    
    // cout << " x1=" << x1 << " k=" << k << " RAreaPeak=" << RAreaPeak << " layer=" << layer << " Ntrue=" << Ntrue << " Npixels=" << Npixels << "  Nmeas/k=" << Nmeas/k << endl;
    return Nmeas/k;    // return in ADC counts
}
Double_t Nonlin_vmeas_func (Double_t *x, Double_t *par) {
    
    // Returns the pulse height V=x1 corrected non-linear correction for fininte number of SiPM pixels.
    
    Double_t k = par[0];            // number of pixels/count
    Double_t RAreaPeak = par[1];            // Ratio of area to peak
    Double_t Npixels = par[2];      // number of pixels in one SiPM array
    Double_t layer = par[3];        // layer number for the correction
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    if (x1 <= 0) return 0;
    
    Npixels *= layer;
    Double_t Ntrue = k*RAreaPeak*x1;
    
    Double_t Nmeas = Npixels*(1-exp(-Ntrue/Npixels));
    Double_t Vmeas = Nmeas/(k*RAreaPeak);
    
    // cout << " x1=" << x1 << " k=" << k << " RAreaPeak=" << RAreaPeak << " layer=" << layer << " Ntrue=" << Ntrue << " Npixels=" << Npixels << "  Nmeas=" << Nmeas << " Vmeas=" << Vmeas << endl;
    return Vmeas;
}
Double_t propo_func (Double_t *x, Double_t *par) {
    
    // Returns funcition y = Ax
    
    Double_t slope = par[0];            // slope
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    // cout << " x1=" << x1 << " slope=" << slope << endl;
    
    if (x1 <= 0) return 0;
    
    return slope*x1;
}
void SiPM_saturation_fits(void)
{
// Pixel non-linearity of the SiPMs plotted as a function of number of pixels firing
//

#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
// gStyle->SetOptStat(kFALSE);
//
    gStyle->SetOptStat(11111111);
    gStyle->SetOptFit(1);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.2);
gStyle->SetPadBottomMargin(0.15);
//
    
    
    Int_t const npar = 4;
    Double_t xmin=0;
    Double_t xmax=4095;
    Double_t ymin=0;
    Double_t ymax=0.3;
    // Double_t k = 0.478;     // number of pixels per count
    Double_t k = 0.6;     // number of pixels per count
    // Double_t RAreaPeak = 12.8;       // ratio of area to peak
    Double_t RAreaPeak = 7.0;       // ratio of area to peak
    Double_t Npixels = 57600;      // number of pixels in one SiPM array
    Double_t layer = 1;        // layer number for the correction
    Double_t k2 = 1;
    
    
    TF1 *Nonlin_vmeas = new TF1 ("Nonlin_vmeas",Nonlin_vmeas_func,xmin,xmax,npar);
    Nonlin_vmeas->SetParameters(k,RAreaPeak,Npixels,layer);
    Nonlin_vmeas->SetParNames("k","RAreaPeak","Npixels","layer");
    
    TF1 *Nonlin_int = new TF1 ("Nonlin_int",Nonlin_int_func,xmin,xmax*RAreaPeak,npar+1);
    k = 0.6;     // number of pixels per count
    // RAreaPeak = 12.8;       // ratio of area to peak
    Npixels = 57600;      // number of pixels in one SiPM array
    layer = 1;        // layer number for the correction
    k2 = 1;
    Nonlin_int->SetParameters(k,RAreaPeak,Npixels,layer,k2);
    Nonlin_int->SetParNames("k","RAreaPeak","Npixels","layer","k2");
    
    const Int_t npts = 100;
    
    for (Int_t jj=1; jj<npts; jj++) {
        Double_t x = xmin + jj*(xmax-xmin)/npts;
        // cout << " jj=" << jj << " k=" << k << " RAreaPeak=" << RAreaPeak << " Npixels=" << Npixels << " layer=" << layer << " x=" << x << " Nonlin_int=" << Nonlin_int->Eval(x) << endl;
    }
    
    
    Double_t fratio1[npts];
    Double_t vtrue[npts];

    for (Int_t j=1; j<npts+1; j++) {
        vtrue[j-1] = j*xmax/npts;
        fratio1[j-1] = vtrue[j-1]>0? Nonlin_vmeas->Eval(vtrue[j-1])/vtrue[j-1] : 0;
        // cout << " j-1=" << j-1 << " layer=" << layer << " vtrue=" << vtrue[j-1] << " fratio1=" << fratio1[j-1] << endl;
    }
    TGraph *gr_fratio1 = new TGraph(npts,vtrue,fratio1);
    
    TCanvas *c1 = new TCanvas("c1","c1 SiPM_saturation_fits",200,10,800,700);
    // c1->Divide(2,2);
    // c1->cd(1);
    // gPad->SetLogy();
  
    layer=1;
    Nonlin_vmeas->SetParameters(k,RAreaPeak,Npixels,layer);
    TF1* Nonlin_vmeas_c1 = Nonlin_vmeas->DrawCopy();
    Nonlin_vmeas_c1->SetTitle("");
    Nonlin_vmeas_c1->GetYaxis()->SetLabelSize(0.05);
    Nonlin_vmeas_c1->GetYaxis()->SetTitleSize(0.07);
    Nonlin_vmeas_c1->GetXaxis()->SetLabelSize(0.05);
    Nonlin_vmeas_c1->GetXaxis()->SetTitleSize(0.07);
    Nonlin_vmeas_c1->GetYaxis()->SetTitleOffset(1.5);
    Nonlin_vmeas_c1->GetYaxis()->SetTitle("Measured Voltage (V)");
    Nonlin_vmeas_c1->GetXaxis()->SetTitle("True Peak(FADC counts)");
    Nonlin_vmeas_c1->GetXaxis()->SetNdivisions(505);
    Nonlin_vmeas_c1->SetLineColor(2);
    Nonlin_vmeas_c1->GetXaxis()->SetRangeUser(xmin,xmax);
    // Nonlin_vmeas_c1->GetYaxis()->SetRangeUser(ymin,ymax);
    // Nonlin_vmeas_c1->Draw();
    
    // c1->cd(4);
    TString string="";
    string.Form("k=%.2f pixel/count",k);
    TLatex *text = new TLatex (1000,3000,string);
    // text->SetNDC();
    text->SetTextSize(0.03);
    text->Draw();
    
    ymin = 0.5;
    ymax = 1.5;
    TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
    
    TCanvas *c2 = new TCanvas("c2","c2 SiPM_saturation_fits",200,10,800,700);
    gr_fratio1->SetTitle("");
    gr_fratio1->GetYaxis()->SetLabelSize(0.05);
    gr_fratio1->GetYaxis()->SetTitleSize(0.07);
    gr_fratio1->GetXaxis()->SetLabelSize(0.05);
    gr_fratio1->GetXaxis()->SetTitleSize(0.07);
    gr_fratio1->GetYaxis()->SetTitleOffset(1.5);
    gr_fratio1->GetYaxis()->SetTitle("Measured / True");
    gr_fratio1->GetXaxis()->SetTitle("True Peak (FADC counts)");
    gr_fratio1->GetXaxis()->SetNdivisions(505);
    gr_fratio1->SetLineColor(2);
    gr_fratio1->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_fratio1->GetYaxis()->SetRangeUser(ymin,ymax);
    
    gr_fratio1->Draw();
    text->DrawLatex(500,1.4,string);
    
    leg->AddEntry(gr_fratio1,"Layer 1","l");
    
    leg->Draw();
    
    // Fit data. From Dhillon Ross and Breanna Crompvoets.  6/25/2019
    
    // Int_t const ndata=9;              // Data from Dhillon and Breanna
    // Double_t filter[ndata]={0.8,0.7,0.63,0.5,0.25,0.10,0.05,0.025,0.0125};                    // Data from Dhillon and Breanna
    // Double_t SiPM_int[ndata]={19440, 18570, 17090, 13500, 7618, 3795, 2052, 1124, 602.8};     // Data from Dhillon and Breanna
    Int_t const ndata=8;              // Data from Dhillon and Breanna
    Double_t filter[ndata]={0.708,0.631,0.5,0.25,0.10,0.05,0.025,0.0125};                    // Log book p 30-31
    Double_t filter_scale[ndata];                    // scaled filter axis
    Double_t filter_err_scale[ndata];                    //  for plotting
    Double_t filter_err[ndata];                    //  for plotting
    Double_t SiPM_int[ndata]={21194, 19410, 16080, 8669, 3947, 2234, 1163, 612};     // Log book p 30-31
    Double_t SiPM_peak[ndata];
    Double_t SiPM_err_peak[ndata];     // Assign a 5% uncertainty due to warmup during measurements
    Double_t SiPM_err_int[ndata];     // Assign a 5% uncertainty due to warmup during measurements
    Double_t DataAreaPeak = 7;
    
    Double_t slope = 1.581*0.937;               // Log p 30-31
    // Double_t slope = 1.14;             // linear fit to 1200
    // Double_t slope = 1.14*1.14;     // linear fit to 800
    // Double_t slope = 1.14*0.88;     // linear fit to 2400
    // convert integral to peak
    for (Int_t j=0; j<ndata; j++) {
        filter_scale[j] = 4000*slope*filter[j];
        filter_err[j] = 0.001;
        filter_err_scale[j] = 1;
        SiPM_peak[j] = SiPM_int[j]/DataAreaPeak;
        SiPM_err_peak[j] = 0.05*SiPM_peak[j];
        SiPM_err_int[j] = 0.05*SiPM_int[j];
    }
    
    TGraph *gr_SiPM_peak = new TGraphErrors(ndata,filter_scale,SiPM_peak,filter_err_scale,SiPM_err_peak);
    
    TCanvas *c3 = new TCanvas("c3","c3 SiPM_saturation_fits",200,10,800,700);
    c3->SetGridx();
    c3->SetGridy();
    
    ymin=0;
    ymax=5000;
    xmin=0;
    xmax=5000;
    gr_SiPM_peak->SetTitle("");
    gr_SiPM_peak->GetYaxis()->SetLabelSize(0.05);
    gr_SiPM_peak->GetYaxis()->SetTitleSize(0.07);
    gr_SiPM_peak->GetXaxis()->SetLabelSize(0.05);
    gr_SiPM_peak->GetXaxis()->SetTitleSize(0.07);
    gr_SiPM_peak->GetYaxis()->SetTitleOffset(1.5);
    string.Form("Integral/%.2f (counts)",DataAreaPeak);
    gr_SiPM_peak->GetYaxis()->SetTitle(string);
    gr_SiPM_peak->GetXaxis()->SetTitle("Scaled Light Intensity (counts)");
    gr_SiPM_peak->GetXaxis()->SetNdivisions(505);
    gr_SiPM_peak->SetLineColor(2);
    gr_SiPM_peak->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_SiPM_peak->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_SiPM_peak->SetMarkerStyle(20);
    gr_SiPM_peak->SetMarkerColor(2);
    gr_SiPM_peak->Draw("Ap");
    
    slope = 1;
    xmin=0;
    xmax=5000;
    TF1 *propo = new TF1 ("propo",propo_func,xmin,xmax,1);
    propo->SetParameter(0,slope);
    propo->SetParName(0,"slope");
    // gr_SiPM_peak->Fit("propo","","",xmin,800);         // linear fit to 800
    // gr_SiPM_peak->Fit("propo","","",xmin,1200);            // linear fit to 1200
    gr_SiPM_peak->Fit("propo","","",xmin,2400);           // linear fit to 2400
    propo->SetLineColor(4);
    propo->DrawCopy("same");
    
    TF1 *fit = gr_SiPM_peak->GetFunction("propo");     // use to scale x-axis to y-axis input
    Double_t a = fit->GetParameter(0);
    Double_t siga = fit->GetParError(0);
    Double_t chi2= fit->GetChisquare();
    Int_t ndf= fit->GetNDF();
    Double_t prob= fit->GetProb();
    
    cout << " a=" << a << " siga=" << siga << " chi2=" << chi2 << endl;
    
    // Nonlin_vmeas->SetParameter(0,4000);
    Nonlin_vmeas->FixParameter(1,RAreaPeak);
    Nonlin_vmeas->FixParameter(2,Npixels);
    Nonlin_vmeas->FixParameter(3,layer);
    gr_SiPM_peak->Fit("Nonlin_vmeas");
    
    TCanvas *c4 = new TCanvas("c4","c4 SiPM_saturation_fits",200,10,800,700);
    c4->SetGridx();
    c4->SetGridy();
    
    // first stab at diode normalization (ignore)
    Double_t filter_pixel[ndata];                    // Log book p 30-31
    Double_t filter_pixel_err[ndata];                    // Log book p 30-31
    
    for (Int_t j=0; j<ndata; j++) {
        filter_pixel[j] = 12400*filter[j];
        filter_pixel_err[j] = 1;
    }
    TGraph *gr_SiPM_int = new TGraphErrors(ndata,filter_pixel,SiPM_int,filter_pixel_err,SiPM_err_int);

    ymin=0;
    ymax=35000;
    xmin=0;
    xmax=15000;
    gr_SiPM_int->SetTitle("");
    gr_SiPM_int->GetYaxis()->SetLabelSize(0.05);
    gr_SiPM_int->GetYaxis()->SetTitleSize(0.07);
    gr_SiPM_int->GetXaxis()->SetLabelSize(0.05);
    gr_SiPM_int->GetXaxis()->SetTitleSize(0.07);
    gr_SiPM_int->GetYaxis()->SetTitleOffset(1.5);
    gr_SiPM_int->GetYaxis()->SetTitle("Integral (counts)");
    gr_SiPM_int->GetXaxis()->SetTitle("Light Intensity (pixels)");
    gr_SiPM_int->GetXaxis()->SetNdivisions(505);
    gr_SiPM_int->SetLineColor(2);
    gr_SiPM_int->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_SiPM_int->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_SiPM_int->SetMarkerStyle(20);
    gr_SiPM_int->SetMarkerColor(2);
    gr_SiPM_int->Draw("Ap");
    Nonlin_int->FixParameter(1,RAreaPeak);
    Nonlin_int->FixParameter(2,Npixels);
    Nonlin_int->FixParameter(3,layer);
    k2 = 1;
    k = 0.3;
    // Nonlin_int->FixParameter(0,k);
    gr_SiPM_int->Fit("Nonlin_int");
    
    TCanvas *c5 = new TCanvas("c5","c5 SiPM_saturation_fits",200,10,800,700);
    c5->SetGridx();
    c5->SetGridy();
    
    TGraph *gr_SiPM_filter = new TGraphErrors(ndata,filter,SiPM_int,filter_err,SiPM_err_int);
    
    ymin=0;
    ymax=35000;
    xmin=0;
    xmax=2;
    gr_SiPM_filter->SetTitle("");
    gr_SiPM_filter->GetYaxis()->SetLabelSize(0.05);
    gr_SiPM_filter->GetYaxis()->SetTitleSize(0.07);
    gr_SiPM_filter->GetXaxis()->SetLabelSize(0.05);
    gr_SiPM_filter->GetXaxis()->SetTitleSize(0.07);
    gr_SiPM_filter->GetYaxis()->SetTitleOffset(1.5);
    gr_SiPM_filter->GetYaxis()->SetTitle("Integral (counts)");
    gr_SiPM_filter->GetXaxis()->SetTitle("Relative filter settings");
    gr_SiPM_filter->GetXaxis()->SetNdivisions(505);
    gr_SiPM_filter->SetLineColor(2);
    gr_SiPM_filter->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_SiPM_filter->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_SiPM_filter->SetMarkerStyle(20);
    gr_SiPM_filter->SetMarkerColor(2);
    gr_SiPM_filter->Draw("Ap");        // linear fit to 800
    slope = 50000;
    propo->SetParameter(0,slope);
    gr_SiPM_filter->Fit("propo","","",xmin,0.07);            // linear fit to 1200
    propo->SetLineColor(4);
    propo->DrawCopy("same");
    
    /*Nonlin_int->FixParameter(1,RAreaPeak);
    Nonlin_int->FixParameter(2,Npixels);
    Nonlin_int->FixParameter(3,layer);
    gr_SiPM_filter->Fit("Nonlin_int");*/
    
    // title.Form("%.0f",rho_dirt*100);
    c4->SaveAs("SiPM_saturation_fits.png");
    // c1->SaveAs("SiPM_saturation_fits.pdf(");
    // c2->SaveAs("SiPM_saturation_fits.pdf");
    c3->SaveAs("SiPM_saturation_fits.pdf(");
    // c4->SaveAs("SiPM_saturation_fits.pdf");
    c5->SaveAs("SiPM_saturation_fits.pdf)");
}