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_lasercalib(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;
    
    // Data from log Book. p 27-28. Take average and std dev of currents (exclude first data set) for measurement.
    
    Int_t const ndata=7;
    Double_t frequency[ndata]={0.1,1,10,100,1000,10000,1000000};
    Double_t frequency_error[ndata]={0.001,0.001,0.001,0.001,0.001,0.001,0.001};
    // Double_t current[ndata]={0.0077,0.0155,0.0902,0.8321,8.5731,96.9250,1000};    // position 1: center of post
    // Double_t current_error[ndata]={0.0004,0.0005,0.0054,0.0596,0.6129,3.0148,1};  //position 1: center of post
    // Double_t current[ndata]={0.0001,0.0075, 0.0415, 0.374, 3.95,1000,1000};    //position 2 (7/3/2019) : side of platform (closest do door) location of A0
    // Double_t current_error[ndata]={0.0001, 0.0001, 0.0011, 0.0120, 0.0354,1,1};  // position 2 (7/3/2019): side of platform (closest do door) location of A0
    Double_t current[ndata]={0.0001,0.0081, 0.0499, 0.4650, 4.7680, 52.3000, 1000};    //position 2 (7/3/2019) : side of platform (closest do door) location of A0
    Double_t current_error[ndata]={0.0001,0.0003, 0.0011, 0.0098, 0.0995, 0.6745,1};  // position 2 (7/3/2019): side of platform (closest do door) location of A0
    Double_t slope = 1;
    
    TGraphErrors *gr_current = new TGraphErrors(ndata,frequency,current,frequency_error,current_error);
    
    TCanvas *c1 = new TCanvas("c1","c1 SiPM_saturation_lasercalib",200,10,800,700);
    c1->SetGridx();
    c1->SetGridy();
    c1->SetLogy();
    c1->SetLogx();
    
    xmin=0.5;
    xmax=20000;
    ymin=0.001;
    ymax=1000;
    gr_current->SetTitle("Position 2");
    gr_current->GetYaxis()->SetLabelSize(0.05);
    gr_current->GetYaxis()->SetTitleSize(0.07);
    gr_current->GetXaxis()->SetLabelSize(0.05);
    gr_current->GetXaxis()->SetTitleSize(0.07);
    gr_current->GetYaxis()->SetTitleOffset(1.5);
    gr_current->GetYaxis()->SetTitle("Current (nA)");
    gr_current->GetXaxis()->SetTitle("Frequency (kHz)");
    gr_current->GetXaxis()->SetNdivisions(505);
    gr_current->SetLineColor(2);
    gr_current->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_current->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_current->SetMarkerStyle(20);
    gr_current->SetMarkerColor(2);
    gr_current->Draw("Ap");
    gr_current->Fit("pol1","","",5,20000);
    // gr_current->Fit("pol1","","",5,50000);
    
    /*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_current->Fit("propo","","",xmin,1200);            // linear fit to 1200
    propo->Draw("same");
    propo->SetLineColor(4);*/
    
    /*TF1 *fit = gr_current->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;*/
    
    
    // title.Form("%.0f",rho_dirt*100);
    c1->SaveAs("SiPM_saturation_lasercalib.png");
    c1->SaveAs("SiPM_saturation_lasercalib.pdf");
}