Double_t Nonlin_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 R = par[0];            // ratio = N(fire) / N(true)
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    if (x1 <= 0) return 0;
    
    Double_t y = (1-exp(-x1))/R;
    return y;    //
}
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_solveR(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 = 1;
    Double_t xmin=0;
    Double_t xmax=0.3;
    Double_t ymin=0;
    Double_t ymax=2;
    // 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
    Double_t k2 = 1;
    
    Double_t R=0.9;     // Ratio = N(fire)/N(true);
    
    TF1 *Nonlin = new TF1 ("Nonlin",Nonlin_func,xmin,xmax,npar);
    Nonlin->SetParameter(0,R);
    Nonlin->SetParName(0,"R");
    
    Double_t slope = 1;
    TF1 *propo = new TF1 ("propo",propo_func,xmin,xmax,1);
    propo->SetParameter(0,slope);
    propo->SetParName(0,"slope");
    
    TCanvas *c1 = new TCanvas("c1","c1 SiPM_saturation_solveR",200,10,800,700);
    c1->SetGridx();
    c1->SetGridy();
    

    ymin=0;
    ymax=35000;
    xmin=0;
    xmax=15000;
    Nonlin->Draw();
    TString text;
    text.Form(" R=%.2f",R);
    Nonlin->SetTitle(text);
    Nonlin->GetYaxis()->SetLabelSize(0.05);
    Nonlin->GetYaxis()->SetTitleSize(0.07);
    Nonlin->GetXaxis()->SetLabelSize(0.05);
    Nonlin->GetXaxis()->SetTitleSize(0.07);
    Nonlin->GetYaxis()->SetTitleOffset(1.5);
    Nonlin->GetYaxis()->SetTitle("y(x,R)");
    Nonlin->GetXaxis()->SetTitle("x");
    Nonlin->GetXaxis()->SetNdivisions(505);
    Nonlin->SetLineColor(2);
    // Nonlin->GetXaxis()->SetRangeUser(xmin,xmax);
    // Nonlin->GetYaxis()->SetRangeUser(ymin,ymax);
    Nonlin->SetLineColor(2);
    Nonlin->SetLineWidth(2);
    
    propo->SetLineColor(4);
    propo->SetLineWidth(2);
    propo->DrawCopy("same");
    
    c1->SaveAs("SiPM_saturation_solveR.png");
    c1->SaveAs("SiPM_saturation_solveR.pdf");
    
}