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 Npixels = 57600;      // number of pixels in one SiPM array
    Double_t layer = 1;        // layer number for the correction
    
    
    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");
    
    const Int_t npts = 100;
    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;
    Double_t filter[ndata]={0.8,0.7,0.63,0.5,0.25,0.10,0.05,0.025,0.0125};
    Double_t SiPM_int[ndata]={19440, 18570, 17090, 13500, 7618, 3795, 2052, 1124, 602.8};
    Double_t DataAreaPeak = 7;
    
    Double_t slope = 1.14;             // linear fit to 1200
    // Double_t slope = 1.14*1.22;     // 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[j] *= 4000*slope;
        SiPM_int[j] /= DataAreaPeak;
    }
    
    TGraph *gr_SiPM_int = new TGraph(ndata,filter,SiPM_int);
    
    TCanvas *c3 = new TCanvas("c3","c3 SiPM_saturation_fits",200,10,800,700);
    c3->SetGridx();
    c3->SetGridy();
    
    ymin=0;
    ymax=5000;
    xmin=0;
    xmax=4000;
    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("Measured");
    gr_SiPM_int->GetXaxis()->SetTitle("Relative Light Intensity");
    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");
    
    slope = 1;
    xmin=0;
    xmax=4000;
    TF1 *propo = new TF1 ("propo",propo_func,xmin,xmax,1);
    propo->SetParameter(0,slope);
    propo->SetParName(0,"slope");
    gr_SiPM_int->SetMarkerStyle(20);
    gr_SiPM_int->SetMarkerColor(2);
    // gr_SiPM_int->Fit("propo","","",xmin,800);         // linear fit to 800
    gr_SiPM_int->Fit("propo","","",xmin,1200);            // linear fit to 1200
    // gr_SiPM_int->Fit("propo","","",xmin,2400);           // linear fit to 2400
    propo->Draw("same");
    propo->SetLineColor(4);
    
    TF1 *fit = gr_SiPM_int->GetFunction("propo");
    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_int->Fit("Nonlin_vmeas");
    
    
    // title.Form("%.0f",rho_dirt*100);
    c3->SaveAs("SiPM_saturation_fits.png");
    c1->SaveAs("SiPM_saturation_fits.pdf(");
    c2->SaveAs("SiPM_saturation_fits.pdf");
    c3->SaveAs("SiPM_saturation_fits.pdf)");
}