Double_t Nonlin_corr_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 Nbar4095 = par[0];     // average number of pixels firing when V=4095
    Double_t Npixels = par[1];      // number of pixels in one SiPM array
    Double_t layer = par[2];        // layer number for the correction
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    //cout << " x1=" << x1 << " par0=" << Nbar4095 << " par1=" << Npixels << " layer=" << layer << endl;
    
    if (x1 <= 0) return 0;
    
    Double_t Nmeas = Nbar4095*(x1/4095);
    if (Nmeas >= Npixels*layer) Nmeas = Npixels*layer -1;   // safeguard against infinity
        
    Double_t Ntrue = -Npixels*layer*log(1-Nmeas/(Npixels*layer));
    Double_t Vtrue = Ntrue * 4095/Nbar4095;
    
    // cout << " x1=" << x1 << " Nmeas=" << Nmeas << " Ntrue=" << Ntrue << " Vtrue=" << Vtrue << endl;
    return Vtrue;
}
void SiPM_saturation_calc(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);
//
    
    
    npar = 3;
    Double_t xmin=0;
    Double_t xmax=4095;
    Double_t ymin=0;
    Double_t ymax=0.3;
    Double_t Nbar4095 = 24000;     // average number of pixels firing when V=4095
    Double_t Npixels = 57600;      // number of pixels in one SiPM array
    Double_t layer = 1;        // layer number for the correction
    
    
    TF1 *Nonlin_corr = new TF1 ("Nonlin_corr",Nonlin_corr_func,xmin,xmax,npar);
    Nonlin_corr->SetParameters(Nbar4095,Npixels,layer);
    Nonlin_corr->SetParNames("Nbar4095","Npixels","layer");
    
    const Int_t npts = 101;
    Double_t fdiff1[npts], fdiff2[npts], fdiff3[npts], fdiff4[npts];
    Double_t vmeas[npts];

    for (Int_t j=0; j<npts; j++) {
        vmeas[j] = j*xmax/(npts-1);
        fdiff1[j] = vmeas[j]>0? Nonlin_corr->Eval(vmeas[j])/vmeas[j] : 0;
        cout << " j=" << j << " layer=" << layer << " vmeas=" << vmeas[j] << " fdiff1=" << fdiff1[j] << endl;
    }
    TGraph *gr_fdiff1 = new TGraph(npts,vmeas,fdiff1);
    
    layer=2;
    Nonlin_corr->SetParameters(Nbar4095,Npixels,layer);
    
    for (Int_t j=0; j<npts; j++) {
        vmeas[j] = j*xmax/(npts-1);
        fdiff2[j] = vmeas[j]>0? Nonlin_corr->Eval(vmeas[j])/vmeas[j] : 0;
        cout << " j=" << j << " layer=" << layer << " vmeas=" << vmeas[j] << " fdiff2=" << fdiff2[j] << endl;
    }
    TGraph *gr_fdiff2 = new TGraph(npts,vmeas,fdiff2);
    
    layer=3;
    Nonlin_corr->SetParameters(Nbar4095,Npixels,layer);
    
    for (Int_t j=0; j<npts; j++) {
        vmeas[j] = j*xmax/(npts-1);
        fdiff3[j] = vmeas[j]>0? Nonlin_corr->Eval(vmeas[j])/vmeas[j] : 0;
        cout << " j=" << j << " layer=" << layer << " vmeas=" << vmeas[j] << " fdiff3=" << fdiff3[j] << endl;
    }
    TGraph *gr_fdiff3 = new TGraph(npts,vmeas,fdiff3);
    
    layer=4;
    Nonlin_corr->SetParameters(Nbar4095,Npixels,layer);
    
    for (Int_t j=0; j<npts; j++) {
        vmeas[j] = j*xmax/(npts-1);
        fdiff4[j] = vmeas[j]>0? Nonlin_corr->Eval(vmeas[j])/vmeas[j] : 0;
        cout << " j=" << j << " layer=" << layer << " vmeas=" << vmeas[j] << " fdiff4=" << fdiff4[j] << endl;
    }
    TGraph *gr_fdiff4 = new TGraph(npts,vmeas,fdiff4);
    
    TString title="";
    // title.Form("#rho_{dirt}=%.2f g/cm^{3}",rho_dirt);
    
    TCanvas *c1 = new TCanvas("c1","c1 SiPM_saturation_calc",200,10,800,700);
    // c1->Divide(2,2);
    // c1->cd(1);
    // gPad->SetLogy();
  
    layer=1;
    Nonlin_corr->SetParameters(Nbar4095,Npixels,layer);
    Nonlin_corr_c1 = Nonlin_corr->DrawCopy();
    Nonlin_corr_c1->SetTitle("");
    Nonlin_corr_c1->GetYaxis()->SetLabelSize(0.05);
    Nonlin_corr_c1->GetYaxis()->SetTitleSize(0.07);
    Nonlin_corr_c1->GetXaxis()->SetLabelSize(0.05);
    Nonlin_corr_c1->GetXaxis()->SetTitleSize(0.07);
    Nonlin_corr_c1->GetYaxis()->SetTitleOffset(1.5);
    Nonlin_corr_c1->GetYaxis()->SetTitle("True Voltage (V)");
    Nonlin_corr_c1->GetXaxis()->SetTitle("Measured Peak(FADC counts)");
    Nonlin_corr_c1->GetXaxis()->SetNdivisions(505);
    Nonlin_corr_c1->SetLineColor(2);
    Nonlin_corr_c1->GetXaxis()->SetRangeUser(xmin,xmax);
    Nonlin_corr_c1->GetYaxis()->SetRangeUser(ymin,ymax);
    
    // c1->cd(4);
    
     /*   TLatex *text = new TLatex (z[0]+50,0.9*ymax,"dump");
    // text->SetNDC();
    text->SetTextSize(0.03);
    text->Draw();
    text->DrawLatex(z[1]+50,0.9*ymax,"Air");
    text->DrawLatex(z[2]+50,0.9*ymax,"Concrete");
    text->DrawLatex(z[3]+50,0.9*ymax,"Dirt");*/
    
    ymin = 1;
    ymax = 1.4;
    TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
    
    TCanvas *c2 = new TCanvas("c2","c2 SiPM_saturation_calc",200,10,800,700);
    gr_fdiff1->SetTitle("");
    gr_fdiff1->GetYaxis()->SetLabelSize(0.05);
    gr_fdiff1->GetYaxis()->SetTitleSize(0.07);
    gr_fdiff1->GetXaxis()->SetLabelSize(0.05);
    gr_fdiff1->GetXaxis()->SetTitleSize(0.07);
    gr_fdiff1->GetYaxis()->SetTitleOffset(1.5);
    gr_fdiff1->GetYaxis()->SetTitle("True / Measured");
    gr_fdiff1->GetXaxis()->SetTitle("Measured Peak (FADC counts)");
    gr_fdiff1->GetXaxis()->SetNdivisions(505);
    gr_fdiff1->SetLineColor(2);
    gr_fdiff1->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_fdiff1->GetYaxis()->SetRangeUser(ymin,ymax);
    
    gr_fdiff1->Draw();
    gr_fdiff2->SetLineColor(4);
    gr_fdiff2->Draw("same");
    gr_fdiff3->SetLineColor(1);
    gr_fdiff3->Draw("same");
    gr_fdiff4->SetLineColor(3);
    gr_fdiff4->Draw("same");
    
    leg->AddEntry(gr_fdiff1,"Layer 1","l");
    leg->AddEntry(gr_fdiff2,"Layer 2","l");
    leg->AddEntry(gr_fdiff3,"Layer 3","l");
    leg->AddEntry(gr_fdiff4,"Layer 4","l");
    
    leg->Draw();
    
    
    // title.Form("%.0f",rho_dirt*100);
    c2->SaveAs("SiPM_saturation_calc.png");
    c2->SaveAs("SiPM_saturation_calc.pdf");
}