#include <inttypes.h>
#include <iostream>
#include <fstream>
#include <utility>
#include <vector>
#include <map>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <cstring>
#include <sstream>

#include "TLatex.h"
#include "TPaveStats.h"
#include "TGraphPainter.h"
#include "TString.h"
#include "TCollection.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMinuit.h"
#include "TKey.h"
#include "TDatime.h"
#include "TAxis.h"
#include "TLine.h"
#include "TTree.h"
#include "TBranch.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TSystem.h"

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 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_fits2(void) {

	//gStyle->Reset();
	//gROOT->SetStyle("Plain");
    gStyle->SetPalette(1,0);
    // gStyle->SetOptStat(kTRUE);
    //gStyle->SetOptStat(kFALSE);
    // gStyle->SetOptFit(kTRUE);
    //gStyle->SetOptFit(kFALSE);
    gStyle->SetOptFit(11111);
    gStyle->SetOptStat(2220);
	/*gStyle->SetStatH(0.10);
	gStyle->SetStatW(0.60);
	gStyle->SetStatX(0.99);
	gStyle->SetStatY(0.99);*/

    gStyle->SetPadRightMargin(0.15);
    gStyle->SetPadLeftMargin(0.2);
    gStyle->SetPadBottomMargin(0.15);
    //gStyle->SetFillColor(0);
    gStyle->SetTitleOffset(1.2, "Y");
  	gStyle->SetTitleSize(0.07,"xyz");
	gStyle->SetTitleSize(0.08,"h");
	gStyle->SetLabelSize(0.07,"xyz");
	gStyle->SetTitleX(0);
	gStyle->SetTitleAlign(13);
	gStyle->SetNdivisions(505,"xy");
    
    
    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
    Double_t k2 = 1;
    
    
    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");
    
    TGraphErrors *gr_SiPM_scan1 = new TGraphErrors("Ch1_Int_537-559.txt","%lg %lg %lg %lg");
    TGraphErrors *gr_SiPM_scan = new TGraphErrors("Ch1_Int_560-585.txt","%lg %lg %lg %lg");
    // use conversion of laser knob to pixels to obtain x-axis values
    
    Double_t q0 = -1288.7;            // offset
    Double_t q1 = 347.011;            // initial slope
    Double_t k0 = 0.264424;           //  exponential normalization
    Double_t k1 = 5.16953;           //  transition
    Double_t kk2 =  7.7799;           //  exponential slope
    Double_t k3 = -1.97623;           //  exponential quadratic term
    
    Int_t const npts1 = gr_SiPM_scan1->GetN();
    cout << " npts1=" << npts1 << endl;
    Double_t xpixel1[npts1];
    Double_t xpixel1_err[npts1];
    Double_t ypixel1[npts1];
    Double_t ypixel1_err[npts1];
    for (Int_t jj=0; jj<npts1; jj++){
        Double_t x = gr_SiPM_scan1->GetX()[jj];
        xpixel1[jj] =  1.34399e+06 -591669*x +84667.3*x*x -3892.38*x*x*x;    // calibration set 1
        // xpixel1[jj] = q0 + q1*x + k0*x*x*exp(kk2*(x-k1)+k3*(x-k1)*(x-k1));  // calibration set 2
        ypixel1[jj] = gr_SiPM_scan1->GetY()[jj];
        xpixel1_err[jj] = gr_SiPM_scan1->GetEX()[jj];
        ypixel1_err[jj] = gr_SiPM_scan1->GetEY()[jj];
        cout << " jj=" << jj << " xpixel1[jj]=" << xpixel1[jj] << " ypixel1[jj]=" << ypixel1[jj] << " ypixel1_err[jj]=" << ypixel1_err[jj] << endl;
    }
    TGraphErrors *gr_SiPM_int1 = new TGraphErrors(npts1,xpixel1,ypixel1,xpixel1_err,ypixel1_err);
    
    
    Int_t const npts = gr_SiPM_scan->GetN();
    cout << " npts=" << npts << endl;
    Double_t xpixel[npts];
    Double_t xpixel_err[npts];
    Double_t ypixel[npts];
    Double_t ypixel_err[npts];
    for (Int_t jj=0; jj<npts; jj++){
        Double_t x = gr_SiPM_scan->GetX()[jj];
        // xpixel[jj] =  1.34399e+06 -591669*x +84667.3*x*x -3892.38*x*x*x;  // calibration set 1
        xpixel[jj] = q0 + q1*x + k0*x*x*exp(kk2*(x-k1)+k3*(x-k1)*(x-k1));  // calibration set 2
        ypixel[jj] = gr_SiPM_scan->GetY()[jj];
        xpixel_err[jj] = gr_SiPM_scan->GetEX()[jj];
        ypixel_err[jj] = gr_SiPM_scan->GetEY()[jj];
        cout << " jj=" << jj << " xpixel[jj]=" << xpixel[jj] << " ypixel[jj]=" << ypixel[jj] << " ypixel_err[jj]=" << ypixel_err[jj] << endl;
    }
    TGraphErrors *gr_SiPM_int = new TGraphErrors(npts,xpixel,ypixel,xpixel_err,ypixel_err);
    
    TCanvas *c0 = new TCanvas("c0","c0 SiPM_saturation_fits",200,10,800,700);
    c0->SetGridx();
    c0->SetGridy();
    
    
    ymin=0;
    ymax=40000;
    xmin=0;
    xmax=30000;
    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->SetMinimum(0);
    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");
    gr_SiPM_int1->SetMarkerStyle(20);
    gr_SiPM_int1->SetMarkerColor(4);
    // gr_SiPM_int1->Draw("samep");
    
    Double_t slope=1;
    TF1 *propo = new TF1 ("propo",propo_func,xmin,xmax,1);
    propo->SetParameter(0,slope);
    propo->SetParName(0,"slope");
    gr_SiPM_int->Fit("propo","","",xmin,2500);            // linear fit to 1200
    propo->SetLineColor(4);
    propo->DrawCopy("same");
    
    // gr_SiPM_int->Fit("pol2");
    Nonlin_int->FixParameter(1,RAreaPeak);
    Nonlin_int->FixParameter(2,Npixels);
    Nonlin_int->FixParameter(3,layer);
    k2 = 1;
    k = 0.478;
    // Nonlin_int->FixParameter(4,k2);
    Nonlin_int->FixParameter(0,k);
    gr_SiPM_int->Fit("Nonlin_int");
    
    /*k = 0.7;
    k2 = Nonlin_int->GetParameter(4);
    TF1 *Nonlin_int2 = (TF1*)Nonlin_int->DrawCopy("same");
    Nonlin_int2->SetParameters(k,RAreaPeak,Npixels,layer,k2);
    Nonlin_int2->SetParNames("k","RAreaPeak","Npixels","layer","k2");
    Nonlin_int2->SetLineColor(4);*/
    
    
    c0->SaveAs("SiPM_saturation_fits2.png");
    c0->SaveAs("SiPM_saturation_fits2.pdf");
}