#include "TCanvas.h"
#include "TGraph.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TAxis.h"
#include "TMath.h"
#include "TAxis.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TPad.h"
#include "TLatex.h"
#include "TApplication.h"

double mean_rms(double* input_array, int n, double& rms)
{
  double mean = 0;
  rms = 0;
  for(int i = 0; i<n; i++)
    {
      mean += input_array[i];
      rms += input_array[i]*input_array[i];
    }
  mean = mean/(double)n;
  rms = rms/(double)n;
  rms = sqrt((rms - mean*mean));
  return mean;
}

//int read_result(char* result_dir, double read_data[6][16], double par[6], double par_err[6], double cond[6])
int read_result(char* result_dir, double read_data[6][16], double cond[6])
{
  char fname[512];

  // Read fit results
  sprintf(fname,"%s/Results.txt",result_dir);
  FILE* fit_results = fopen(fname,"r");
  if(!fit_results)
    {
      cout<<"results file not found in "<<result_dir<<endl;
      return 0;
    }

  int nl = 0;
  char line[512];
  double tmp_read[6];
  //  double read_data[6][16];

  int nc = 16;
  int np = 6;

  fgets(line,512,fit_results);
  while(fgets(line,512,fit_results)||nl<nc)
    {
      if(sscanf(line, "%*d%lf%*lf%lf%*lf%lf%*lf%lf%*lf%lf%*lf%lf%*lf",
		tmp_read,tmp_read+1,tmp_read+2,tmp_read+3,tmp_read+4,tmp_read+5)==6){
	for(int i = 0; i<np; i++) read_data[i][nl] = tmp_read[i];
	nl++;
      }
    }
  if(nl!=nc) return 0;

  // Read average voltage, tmperature and current
  sprintf(fname,"%s/Average.txt",result_dir);
  FILE* condition = fopen(fname,"r");
  if(!condition)
    {
      cout<<"condition file not found in "<<result_dir<<endl;
      return 0;
    }
  
  double tmpv;
  double norm_voltage = 0;
  double set_voltage = 0;
  double temperature = 0;
  double current = 0;
  double gate = 0;

  fgets(line,512,condition);
  if(sscanf(line, "Average Voltage: %lf", &tmpv)) set_voltage = tmpv;

  fgets(line,512,condition);
  if(sscanf(line, "Average Current: %lf", &tmpv)) current = tmpv;

  fgets(line,512,condition);
  if(sscanf(line, "Average Temperature: %lf", &tmpv)) temperature = tmpv;
  
  fgets(line,512,condition);
  if(sscanf(line, "Bias at 25C: %lf", &tmpv)) norm_voltage = tmpv;

  // Read ADC gate length
  sprintf(fname,"%s/Settings.txt",result_dir);
  FILE* setting = fopen(fname,"r");
  if(!setting)
    {
      cout<<"setting file not found in "<<result_dir<<endl;
      return 0;
    }
  while(fgets(line,512,setting))
    {
      if(sscanf(line, "Gate Width = %lf", &tmpv)){
	gate = tmpv;
	break;
      }
    }

  cond[0] = norm_voltage;
  cond[1] = temperature;
  cond[2] = set_voltage;
  cond[3] = current;
  cond[4] = gate;
  
  //  for(int i = 0; i<np; i++) par[i] = mean_rms(read_data[i],nc,par_err[i]);
  return 1;
}

void extract()
{
  gStyle->SetOptStat(1);

  char data_dir[512];
  int input_offset = 2;
  sprintf(data_dir,"./results/%s",gApplication->Argv(input_offset));
  int correction_flag = atoi(gApplication->Argv(1+input_offset));
  int save_flag = atoi(gApplication->Argv(2+input_offset));

  cout<<"Reading data in "<<data_dir<<endl;
  cout<<"Correct temperature: "<<correction_flag<<endl;
  cout<<"Save result: "<<save_flag<<endl;

  int nf = 5;
  int np = 6;
  int nc = 16;
  char filter[5][20] = {"DarkRun","Filter=1%","Filter=2%","Filter=4%","Filter=6%"};

  // Parameters: npe, gain, ped, xt, sig, sig_p
  double data[16][6][5]; // channels // parameters // filter
  double par[6][5]; // parameters // filter
  double par_err[6][5]; // parameter errors // filter
  double condition[6][5]; // vol, temp, set_v, curr, gate // filter
  double uniformity[4][16]; // variables: dark noise, PDE, gain, x-talk // channels 
  double uni_mean[4], uni_rms[4];

  // Temparory variables
  double tmp_data[6][16],tmp_par[16],tmp_par_e[6],tmp_cond[6];
  char result_dir[128];

  // read data and condition
  for(int i = 0; i<nf; i++)
    {
      sprintf(result_dir,"%s/%s",data_dir,filter[i]);
      if(read_result(result_dir,tmp_data,tmp_cond))
	for(int j = 0; j<np; j++)
	  {
	    condition[j][i] = tmp_cond[j];
	    for(int k = 0; k<nc; k++) 
	      {
		// traslate gain from ADC channel to actual value
		//		if(j == 1) data[k][j][i] = tmp_data[j][k]*0.125; // 0.125x10^5 // 0.936
		if(j == 1) data[k][j][i] = tmp_data[j][k]*0.0936; // 0.125x10^5 // 0.0936
		else data[k][j][i] = tmp_data[j][k];
	      }
	  }
    }

  // Apply Temperature Corrections
  // Correction constants
  double v_t = 0.0564; // bias voltage with temperature: linear : Fig.[6]
  double dr_t = 0.0767; // dark rate with bias: exp : Fig[10]
  double dr_v = 1.155; // dark rate with temperature: exp : Fig.[11]
  double g_v = 5.5; // gain (x10^5) with bias: linear : Fig.[5]
  double vol,temp,set_v,t_off,v_off;

  double p_v = 0.913; // pde vs bias voltage
  double x_v = 0.352; // x-talk vs bias voltage

  if(correction_flag)
    for(int i = 0; i<nc; i++)
      for(int j = 0; j<nf; j++)
	{
	  vol = condition[0][j];
	  temp = condition[1][j];
	  set_v = condition[2][j];
	  t_off = temp-25;
	  v_off = set_v - (vol + v_t*t_off);
	  
	  // dark_rate and pde
	  // pure dark
	  if(j == 0) data[i][0][j] *= exp(-dr_v*v_off)*exp(-dr_t*t_off);
	  else
	    {
	      // correct for dark
	      data[i][0][j] += data[i][0][0]*(1.-exp(dr_v*v_off)*exp(dr_t*t_off));
	      // correct for pde
	      data[i][0][j] = data[i][0][0] + (data[i][0][j] - data[i][0][0])/(1+p_v*v_off);
	    }
	  // gain
	  data[i][1][j] += -g_v*v_off;
	  // xtalk: no way to correct it
	}
  
  // Plot averaged parameters
  
  for(int i = 0; i<nf; i++)
    for(int j = 0; j<np; j++)
      {
	for(int k = 0; k<nc; k++) tmp_data[j][k] = data[k][j][i];
	par[j][i] = mean_rms(tmp_data[j],nc,par_err[j][i]);
      }

  TCanvas* c1 = new TCanvas("c1","",800,960);
  c1->Divide(2,3);
  
  // light intensity: 12.78 photons/mm^2, for each tile 9 mm^2 => 115 photon/tile,
  // percentage
  double l_input[5] = {0., 0.0160,0.0269,0.0528,0.0778};
  for(int i = 0; i<nf; i++) l_input[i] = l_input[i]*115.;
  double l_input_e[5] = {0., 0., 0., 0., 0.};

  TGraphErrors* gr[6];
  Color_t col[6] = {kCyan+2,kRed-3,kAzure+3,kOrange+7,kMagenta-2,kTeal+3};
  char title[6][128] = {"Pixels fired per tile", "Charge of a pixel", "Pedestal", "Cross-Talk + Afterpulsing", "Sigma of single pixel signal", "Sigma of Pedestal"};
  char ytitle[6][128] = {"Number of pixels", "ADC channels", "ADC channels", "X-Talk+AP", "ADC channels", "ADC channels"};
  char grtitle[128];

  for(int i = 0; i<np; i++)
    {
      gr[i] =  new TGraphErrors(nf,l_input,par[i],l_input_e,par_err[i]);
      gr[i]->SetMarkerStyle(20+i);
      gr[i]->SetMarkerColor(col[i]);
      gr[i]->SetLineColor(col[i]);
      gr[i]->GetXaxis()->SetLimits(-1,10);
      sprintf(grtitle,"%s; Number of Photons; %s",title[i],ytitle[i]);
      gr[i]->SetTitle(grtitle);
      c1->cd(i+1);
      gPad->SetGridx(1);
      gPad->SetGridy(1);
      gr[i]->Draw("AP");
      if(i == 0) 
	{
	  gr[i]->Fit("pol1","Q");
	  gr[i]->GetFunction("pol1")->SetLineColor(col[i]);
	  gr[i]->GetFunction("pol1")->SetLineWidth(2);
	}
      else 
	{
	  gr[i]->Fit("pol0","Q");
	  gr[i]->GetFunction("pol0")->SetLineColor(col[i]);
	  gr[i]->GetFunction("pol0")->SetLineWidth(2);
	}
    }

  // Fit Uniformity Variables
  
  TCanvas* c2 = new TCanvas("c2","",1200,800);
  c2->Divide(4,4);

  // get parameters for each channels

  TGraph* tmp_gr[16];
  double tmp_rms;

  for(int i = 0; i<nc; i++)
    {  
      c2->cd(i+1);
 
      // dark noise and PDE
      tmp_gr[i] = new TGraph(nf, l_input, data[i][0]);
      sprintf(grtitle,"Ch#%02d; Number of photons; Number of Pixels",i);
      tmp_gr[i]->SetTitle(grtitle);
      tmp_gr[i]->GetXaxis()->SetLimits(-1,10);
      tmp_gr[i]->SetMarkerStyle(4);
      tmp_gr[i]->SetMarkerColor(2);
      tmp_gr[i]->Fit("pol1","Q");
      tmp_gr[i]->GetFunction("pol1")->SetLineWidth(1);

      // Chage hits to rate by dividing gate
      uniformity[0][i] = tmp_gr[i]->GetFunction("pol1")->GetParameter(0)*1000./condition[4][0];
      uniformity[1][i] = tmp_gr[i]->GetFunction("pol1")->GetParameter(1); // PDE
      uniformity[2][i] = mean_rms(data[i][1],nf,tmp_rms); // gain
      uniformity[3][i] = mean_rms(data[i][3],nf,tmp_rms); // x-talk
      tmp_gr[i]->Draw("AP");
    }
  
  for(int i = 0; i<4; i++) uni_mean[i] = mean_rms(uniformity[i],nc,uni_rms[i]);

  // Plot Uniformities
  TCanvas* c3 = new TCanvas("c3","",800,600);
  c3->Divide(2,2);

  TGraph* bar_gr[4];
  double channel[16];
  double norm_uni[16];

  for(int i = 0; i<nc; i++) channel[i] = (double)i;

  Color_t col_bar[4] = {kRed-4,kAzure+5,kTeal-7,kOrange+7};
  char bar_title[4][256] = {"Dark Rate","PDE","Gain","Cross-Talk"};
  char bar_unit[4][256] = {"MHz","","#times 10^{5}",""};
  TLatex ltext;
  ltext.SetTextSize(0.08);

  for(int i = 0; i<4; i++)
    {
      c3->cd(i+1);
      gPad->SetGridy(1);

      //      for(int j = 0; j<nc; j++) norm_uni[j] = uniformity[i][j]/uni_mean[i];
      for(int j = 0; j<nc; j++) norm_uni[j] = uniformity[i][j]/uni_mean[i]-1.;

      bar_gr[i] = new TGraph(nc, channel, norm_uni);

      sprintf(grtitle,"%s Uniformity; Channel",bar_title[i]);
      bar_gr[i]->SetTitle(grtitle);
      bar_gr[i]->SetFillColor(col_bar[i]);
      bar_gr[i]->GetXaxis()->SetLimits(-1,nc);
      bar_gr[i]->Draw("AB");
      if(i == 1||i ==2)
	{
	  bar_gr[i]->SetMinimum(-0.15);
	  bar_gr[i]->SetMaximum(0.15);
	}
      else
	{
	  bar_gr[i]->SetMinimum(-0.4);
	  bar_gr[i]->SetMaximum(0.4);
	}
      sprintf(grtitle,"%.3f #pm %.1f%% %s",uni_mean[i],uni_rms[i]*100/uni_mean[i],bar_unit[i]);
      if(i == 1||i ==2)
	{
	  ltext.DrawLatex(nc*.2, 0.12, grtitle);
	}
      else
	ltext.DrawLatex(nc*.2, 0.32, grtitle);
      
      c3->Update();
    }
    
  // Save data
  char tale[2][20] = {"_nc",""};
  if(save_flag)
    {
      char fname[512];
      sprintf(fname,"%s/Parameters%s.png",data_dir,tale[correction_flag]);
      c1->Print(fname);

      sprintf(fname,"%s/PDE_Fit%s.png",data_dir,tale[correction_flag]);
      c2->Print(fname);

      sprintf(fname,"%s/Uniformity%s.png",data_dir,tale[correction_flag]);
      c3->Print(fname);

      sprintf(fname,"%s/Summary%s.txt",data_dir,tale[correction_flag]);
      FILE* summary = fopen(fname,"w");

      fprintf(summary,"PhotonInput\t npe\t npe_e\t gain(x10^5)\t gain_e\t pedestal\t ped_e\t xtalk\t xtalk_e\t sigma\t sig_e\t sigma_p\t sig_p_e\n");
      for(int i = 0; i<nf; i++)
	{
	  fprintf(summary, "%.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\t %.6f\n",l_input[i],par[0][i],par_err[0][i],par[1][i],par_err[1][i],par[2][i],par_err[2][i],par[3][i],par_err[3][i],par[4][i],par_err[4][i],par[5][i],par_err[5][i]);
	}
      fprintf(summary,"\nChannel#\tDarkRate (MHz)\tPDE\tGain(x10^5)\tCross-Talk\n");
      for(int i = 0; i<nc; i++)
	{
	  fprintf(summary,"%d\t%.6f\t%.6f\t%.6f\t%.6f\n",i,uniformity[0][i],uniformity[1][i],uniformity[2][i],uniformity[3][i]);
	}
      fprintf(summary,"\nmean\t%.6f\t%.6f\t%.6f\t%.6f\n",uni_mean[0],uni_mean[1],uni_mean[2],uni_mean[3]);
      fprintf(summary,"rms\t%.6f\t%.6f\t%.6f\t%.6f\n",uni_rms[0],uni_rms[1],uni_rms[2],uni_rms[3]);

      char dir_name[50];
      sscanf(data_dir,"./results/%s",dir_name);

      sprintf(fname,"./results/darkrate.txt");
      FILE* darkrate_out = fopen(fname,"a");
      sprintf(fname,"./results/pde.txt");
      FILE* pde_out = fopen(fname,"a");
      sprintf(fname,"./results/gain.txt");
      FILE* gain_out = fopen(fname,"a");
      sprintf(fname,"./results/xtalk.txt");
      FILE* xtalk_out = fopen(fname,"a");

      //summary
      sprintf(fname,"./results/summary.txt");
      FILE* summary_out = fopen(fname,"a");
      fprintf(summary_out,"%s",dir_name);
      fprintf(summary_out,"\t%.6f\t%.6f",condition[0][0],condition[4][0]);
      fprintf(summary_out,"\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\t%.6f\n",uni_mean[0],uni_rms[0]*100./uni_mean[0],uni_mean[1],uni_rms[1]*100./uni_mean[1],uni_mean[2],uni_rms[2]*100./uni_mean[2],uni_mean[3],uni_rms[3]*100./uni_mean[3]);

      fprintf(darkrate_out,"%s",dir_name);
      fprintf(pde_out,"%s",dir_name);
      fprintf(gain_out,"%s",dir_name);
      fprintf(xtalk_out,"%s",dir_name);

      for(int i = 0; i<nc; i++)
	{
	  fprintf(darkrate_out,"\t%.6f",uniformity[0][i]);
	  fprintf(pde_out,"\t%.6f",uniformity[1][i]);
	  fprintf(gain_out,"\t%.6f",uniformity[2][i]);
	  fprintf(xtalk_out,"\t%.6f",uniformity[3][i]);
	}

      fprintf(darkrate_out,"\n");
      fprintf(pde_out,"\n");
      fprintf(gain_out,"\n");
      fprintf(xtalk_out,"\n");

    }
}