#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 knob_calib_func (Double_t *x, Double_t *par) {
    
    // function to calibrate the photon flux vs laser knob setting.
    
    Double_t q0 = par[0];            // offset
    Double_t q1 = par[1];            // initial slope
    Double_t k0 = par[2];           //  exponential normalization
    Double_t k1 = par[3];           //  transition
    Double_t k2 = par[4];           //  exponential slope
    Double_t k3 = par[5];           //  exponential quadratic term
    Double_t x1 = x[0] ;            // input pulse height in counts
    
    if (x1 <= 0) return 0;
    
    Double_t f = q0 + q1*x1 + k0*x1*x1*exp(k2*(x1-k1)+k3*(x1-k1)*(x1-k1));
    
    // cout << " x1=" << x1 << " k=" << k << " RAreaPeak=" << RAreaPeak << " layer=" << layer << " Ntrue=" << Ntrue << " Npixels=" << Npixels << "  Nmeas/k=" << Nmeas/k << endl;
    return f;    // return in ADC counts
}

void SiPM_saturation_knobcalib() {

	//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 = 6;
    Double_t xmin=4;
    Double_t xmax=8;
    Double_t ymin=0;
    Double_t ymax=2000;
    Double_t q0=5000;     // offset
    Double_t q1=-400;       // slope at zero
    Double_t k0=0.02;      // exponential norm
    Double_t k1=5;        //  transition
    Double_t k2 = 7;        // exponential slope
    Double_t k3 = 0;        // exponential slope
    
    
    TF1 *knob_calib = new TF1 ("knob_calib",knob_calib_func,xmin,xmax,npar);
    knob_calib->SetParameters(q0,q1,k0,k1,k2,k3);
    knob_calib->SetParNames("q0","q1","k0","k1","k2","k3");
    

  TCanvas *canvasall = new TCanvas("canvasall","canvasall",1000,800);
  canvasall->Divide(3,2,0.001,0.001);
    
  /*canvasall->cd(1);
  TGraphErrors *gr_f10kHz = new TGraphErrors("photodiode_calib_10kHz.txt","%lg %lg %lg");
    gr_f10kHz->SetTitle("Set 1: 10kHz;laser knob setting;diode current (pA)");
  gr_f10kHz->SetMarkerStyle(20);
  gr_f10kHz->SetMinimum(0);
  gr_f10kHz->Draw("ap");*/

  canvasall->cd(2);
  TGraphErrors *gr_f100kHz = new TGraphErrors("photodiode_calib2_100kHz.txt","%lg %lg %lg");
    gr_f100kHz->SetTitle("Set 1: 100kHz;laser knob setting;diode current (pA)");
  gr_f100kHz->SetMarkerStyle(20);
  gr_f100kHz->SetMinimum(0);
  gr_f100kHz->Draw("ap");
    
    canvasall->cd(3);
    TGraphErrors *gr_f50kHz = new TGraphErrors("photodiode_calib2_50kHz.txt","%lg %lg %lg");
    gr_f50kHz->SetTitle("Set 1: 50kHz;laser knob setting;diode current (pA)");
    gr_f50kHz->SetMarkerStyle(20);
    gr_f50kHz->SetMinimum(0);
    gr_f50kHz->Draw("ap");
    
    /*canvasall->cd(4);
    TGraphErrors *gr2_f100kHz = new TGraphErrors("photodiode_calib2_100kHz.txt","%lg %lg %lg");
    gr2_f100kHz->SetTitle("Set 2: 100kHz;laser knob setting;diode current (pA)");
    gr2_f100kHz->SetMarkerStyle(20);
    gr2_f100kHz->SetMinimum(0);
    gr2_f100kHz->Draw("ap");
    
    canvasall->cd(5);
    TGraphErrors *gr2_f50kHz = new TGraphErrors("photodiode_calib2_50kHz.txt","%lg %lg %lg");
    gr2_f50kHz->SetTitle("Set 2: 50kHz;laser knob setting;diode current (pA)");
    gr2_f50kHz->SetMarkerStyle(20);
    gr2_f50kHz->SetMinimum(0);
    gr2_f50kHz->Draw("ap");*/
    
  // get data and scale to flux
    // Scale to nominal flux for laser knob setting of 7.0: 54523 photons/pulse/100 mm2 -> 78513 photons/pulse/144 mm2
    
    
    // Double_t offset10=2.74;
    // Double_t offset50=2.74;
    // Double_t offset100=2.74;
    Double_t offset50=2.8;
    Double_t offset100=2.8;
    // Double_t SiQE=0.2;      // initial guess
    Double_t SiQE=0.356;      // From JLab data base: PDE=0.246 (V=0.9) x 1.45 = 0.356
    Double_t photons_per_pulse = 54066 * 1.44;     // scale by area (SiPM/photodiode)
    Double_t pixels_per_pulse = SiQE*photons_per_pulse;
    // Double_t scale10=pixels_per_pulse/(gr_f10kHz->GetY()[0]-offset10);
    Double_t scale50=pixels_per_pulse/(gr_f50kHz->GetY()[0]-offset50);
    Double_t scale100=pixels_per_pulse/(gr_f100kHz->GetY()[0]-offset100);
    
    /*Int_t const npts10 = gr_f10kHz->GetN();
    cout << " npts10=" << npts10 << endl;
    Double_t x10kHz[npts10];
    Double_t x10kHz_err[npts10];
    Double_t y10kHz[npts10];
    Double_t y10kHz_err[npts10];
    Double_t y10kHz_scale[npts10];
    Double_t y10kHz_err_scale[npts10];
    for (Int_t jj=0; jj<npts10; jj++){
        x10kHz[jj] = gr_f10kHz->GetX()[jj];
        y10kHz[jj] = gr_f10kHz->GetY()[jj]-offset10;
        y10kHz_scale[jj] = scale10*y10kHz[jj];
        x10kHz_err[jj] = 0;
        y10kHz_err[jj] = gr_f10kHz->GetEY()[jj];
        y10kHz_err_scale[jj] = scale10*y10kHz_err[jj];
        cout << " jj=" << jj << " x10kHz[jj]=" << x10kHz[jj] << " y10kHz[jj]=" << y10kHz_scale[jj] << " y10kHz_err[jj]=" << y10kHz_err_scale[jj] << endl;
    }*/
    
    Int_t const npts50 = gr_f50kHz->GetN();
    cout << " npts50=" << npts50 << endl;
    Double_t x50kHz[npts50];
    Double_t x50kHz_err[npts50];
    Double_t y50kHz[npts50];
    Double_t y50kHz_err[npts50];
    Double_t y50kHz_scale[npts50];
    Double_t y50kHz_err_scale[npts50];
    for (Int_t jj=0; jj<npts50; jj++){
        x50kHz[jj] = gr_f50kHz->GetX()[jj];
        y50kHz[jj] = gr_f50kHz->GetY()[jj]-offset50;
        y50kHz_scale[jj] = scale50*y50kHz[jj];
        x50kHz_err[jj] = 0;
        y50kHz_err[jj] = gr_f50kHz->GetEY()[jj];
        y50kHz_err_scale[jj] = scale50*y50kHz_err[jj];
        cout << " jj=" << jj << " x50kHz[jj]=" << x50kHz[jj] << " y50kHz[jj]=" << y50kHz_scale[jj] << " y50kHz_err[jj]=" << y50kHz_err_scale[jj] << endl;
    }
    
    Int_t const npts100 = gr_f100kHz->GetN();
    cout << " npts100=" << npts100 << endl;
    Double_t x100kHz[npts100];
    Double_t x100kHz_err[npts100];
    Double_t y100kHz[npts100];
    Double_t y100kHz_err[npts100];
    Double_t y100kHz_scale[npts100];
    Double_t y100kHz_err_scale[npts100];
    for (Int_t jj=0; jj<npts100; jj++){
        x100kHz[jj] = gr_f100kHz->GetX()[jj];
        y100kHz[jj] = gr_f100kHz->GetY()[jj]-offset100;
        y100kHz_scale[jj] = scale100*y100kHz[jj];
        x100kHz_err[jj] = 0;
        y100kHz_err[jj] = gr_f100kHz->GetEY()[jj];
        y100kHz_err_scale[jj] = scale100*y100kHz_err[jj];
        cout << " jj=" << jj << " x100kHz[jj]=" << x100kHz[jj] << " y100kHz[jj]=" << y100kHz_scale[jj] << " y100kHz_err[jj]=" << y100kHz_err_scale[jj] << endl;
    }
    
    // Create new scaled TGraphErrors and plot
    
    // TGraphErrors *gr_f10kHz_scaled = new TGraphErrors(npts10,x10kHz,y10kHz_scale,x10kHz_err,y10kHz_err_scale);
    TGraphErrors *gr_f50kHz_scaled = new TGraphErrors(npts50,x50kHz,y50kHz_scale,x50kHz_err,y50kHz_err_scale);
    TGraphErrors *gr_f100kHz_scaled = new TGraphErrors(npts100,x100kHz,y100kHz_scale,x100kHz_err,y100kHz_err_scale);
    
    TCanvas *c0 = new TCanvas("c0","c0 SiPM_saturation_fits",200,10,800,700);
    c0->SetGridx();
    c0->SetGridy();
    
    ymin=0;
    ymax=60000;
    xmin=0;
    xmax=2;
    gr_f100kHz_scaled->SetTitle("");
    gr_f100kHz_scaled->GetYaxis()->SetLabelSize(0.05);
    gr_f100kHz_scaled->GetYaxis()->SetTitleSize(0.07);
    gr_f100kHz_scaled->GetXaxis()->SetLabelSize(0.05);
    gr_f100kHz_scaled->GetXaxis()->SetTitleSize(0.07);
    gr_f100kHz_scaled->GetYaxis()->SetTitleOffset(1.5);
    gr_f100kHz_scaled->GetYaxis()->SetTitle("Pixels per pulse (no saturation)");
    gr_f100kHz_scaled->GetXaxis()->SetTitle("Laser knob setting");
    gr_f100kHz_scaled->GetXaxis()->SetNdivisions(505);
    gr_f100kHz_scaled->SetMinimum(0);
    gr_f100kHz_scaled->SetLineColor(2);
    // gr_f100kHz_scaled->GetXaxis()->SetRangeUser(xmin,xmax);
    gr_f100kHz_scaled->GetYaxis()->SetRangeUser(ymin,ymax);
    gr_f100kHz_scaled->SetMarkerStyle(24);
    gr_f100kHz_scaled->SetMarkerColor(2);
    gr_f100kHz_scaled->Draw("Ap");
    // gr_f100kHz_scaled->Fit("pol3");
    /*gr_f10kHz_scaled->SetMarkerStyle(24);
    gr_f10kHz_scaled->SetMarkerColor(4);
    gr_f10kHz_scaled->Draw("samep");*/
    gr_f50kHz_scaled->SetMarkerStyle(25);
    gr_f50kHz_scaled->SetMarkerColor(2);
    gr_f50kHz_scaled->Fit("knob_calib");
    gr_f50kHz_scaled->Draw("samep");
    
    TF1 *fit = gr_f50kHz_scaled->GetFunction("knob_calib");
    q0 = fit->GetParameter(0);
    Double_t sigq0 = fit->GetParError(0);
    q1 = fit->GetParameter(1);
    Double_t sigq1 = fit->GetParError(1);
    k0 = fit->GetParameter(2);
    Double_t sigk0 = fit->GetParError(2);
    k1 = fit->GetParameter(3);
    Double_t sigk1 = fit->GetParError(3);
    k2 = fit->GetParameter(4);
    Double_t sigk2 = fit->GetParError(4);
    k3 = fit->GetParameter(5);
    Double_t sigk3 = fit->GetParError(5);
    Double_t chi2= fit->GetChisquare();
    Int_t ndf= fit->GetNDF();
    Double_t prob= fit->GetProb();
    
    cout << " par=" << q0 << " " << q1 << " " << k0 << " " << k1 << " " << k2 << " " << k3 << " chi2=" << chi2 << " ndf=" << ndf << endl;
    
    TLegend *leg = new TLegend(0.2,0.7,0.4,0.9);
    // leg->AddEntry(gr_f10kHz_scaled,"10 KHz","p");
    leg->AddEntry(gr_f50kHz_scaled,"50 KHz","p");
    leg->AddEntry(gr_f100kHz_scaled,"100 KHz","p");
    leg->Draw();
    
    c0->SaveAs("SiPM_saturation_knobcalib.png");
    c0->SaveAs("SiPM_saturation_knobcalib.pdf");
    
    
}