// to analyse the backplane pulsing data per module
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)

#include "TApplication.h"
#include "sls_detector_calibration/jungfrauCommonHeader.h"
#include "sls_detector_calibration/jungfrauCommonFunctions.h"

#include "sls_detector_calibration/jungfrauFile.C"
#include "sls_detector_calibration/jungfrauPedestal.C"

#include "TGraphErrors.h"
#include "TF1.h"
#include "TStyle.h"
#include "TPaveStats.h"
#include "TFile.h"
#include "TLegend.h"
#include "TPaveText.h"
#include "TCanvas.h"
#include "TRootCanvas.h"
#include "TSystem.h"
#include "TF2.h"

#include <sys/stat.h>
#include <sstream>

//#define NB_ENABLE 1
//void nonblock(int state);
//#define NB_DISABLE 0

//TApplication* rootapp;
//TCanvas *A2;
//TCanvas *A3;
//TCanvas *A4;
//TCanvas *A5;
//TCanvas *A6;

TGraphErrors *grap_g0;
TGraphErrors *grap_g1;


TF1 *fit_g0=0;
TF1 *fit_g1=0;

TGraphErrors *norm_g0=0;
TGraphErrors *norm_g1=0;

TF1 *flat_g0;
TF1 *flat_g1;

TF1 *lin_g0_p1pc;
TF1 *lin_g0_p05pc;
TF1 *lin_g0_p02pc;
TF1 *lin_g0_m1pc;
TF1 *lin_g0_m05pc;
TF1 *lin_g0_m02pc;

TF1 *lin_g1_p02pc;
TF1 *lin_g1_p01pc;
TF1 *lin_g1_m02pc;
TF1 *lin_g1_m01pc;

//TPaveStats *st0;


//void PlotCanvas(void);

double checkRangeMaxForAmplifierPlateau(double range_max) {

  // check that the range maximum is no more than 6.4 V
  // to avoid non-linearity coming from amplifier plateau
       if (range_max > 6400) {
         return 6400;
       } else {
         return range_max;
         }
}
 



bool isHGX=false;
int main(int argc, char* argv[]) {

 //nonblock(NB_ENABLE);
 cout <<"opening the rootapp" <<endl;
 TApplication rootapp("example",&argc, argv);
  
  jungfrauStyle();
  //gROOT->SetBatch(1);
  gStyle->SetOptFit(11);

  /*
  if (argc != 3) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify data location" << endl;
    cout << "arg 3: specify column (x)" << endl;
    cout << "arg 4: specify row (y)" << endl;
    exit(1);
  }
  */ //uncomment for SR

  if (argc != 4) {  
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify pixel x position" << endl;
    cout << "arg 3: specify pixel y position" << endl;
    exit(1);
  } //uncomment for VH 210906

  string module_str = argv[1];

  string str2 =("HG0G1G2");
  string str3 =("HGOG1G2"); //filename creator had this bug 
   
  string C = argv[2]; 
  string R = argv[3];
  
  
  int column;
  int row;


  std::stringstream(C) >> column;
  std::stringstream(R) >> row;

  int pixel = column+row*1024;

  char data_loc[256];
  sprintf(data_loc,"/mnt/sls_det_storage/jungfrau_calib/jungfrau_ana_sophie/M%s_CalibAna/", module_str.c_str()); 
  cout << data_loc << endl;
  
  std::string folder_path(data_loc);
  if (folder_path.find(str2) != string::npos) isHGX=true;  
  if (folder_path.find(str3) != string::npos) isHGX=true; 


  // cout << data_loc.find(str2)<<" "  << string::npos << " " << str2 << " " << data_loc <<endl;


  if (isHGX) {
    cout << " HG0->HG1->HG2 sequence - dynamicHG0" <<endl;
   // plotfolder_str="BackplanePulsing_HG0G1G2"; 
  }
  else {
    cout << " G0->G1->G2 sequence - dynamicG0" <<endl;
   // plotfolder_str="BackplanePulsing"; 
  }
 
 
  std::vector<double> G0_pixel(245, 0);
  std::vector<double> G0_err_pixel(245, 0);
  std::vector<double> G1_pixel(245, 0);
  std::vector<double> G1_err_pixel(245, 0);  

  
  //char savename[128]; //uncomment for SR
  //char filename[128]; //uncomment for SR
  char filename[256]; //uncomment for VH 210902

  // create necessary directories with permissions drwxrwxr-x
  // data/Mxxx
  //sprintf(savename,"data/M%s", module_str.c_str());
  //mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx
  //sprintf(savename,"plots/M%s", module_str.c_str());
  //mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/BackplanePulsing
  //sprintf(savename,"plots/M%s/%s", module_str.c_str(), plotfolder_str.c_str());
  //mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // /mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/Mxxx_CalibAna
  //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna", module_str.c_str()); //uncomment for SR
  //sprintf(savename,"%s", anadata_loc.c_str()); //uncomment for VH 210906
  //mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  double xs[245];
  for (int i = 0; i < 139; i++) {
    xs[i] = 10+i*5.;
  }
  for (int i = 0; i < 106; i++) {
    xs[i+139] = 700+((i+1)*50);
  }

TCanvas  *A2 = new TCanvas("A2","Plot scan",150,10,800,400);
TCanvas  *A3 = new TCanvas("A3","Plot G0 fit",150,10,800,400);
TCanvas  *A4 = new TCanvas("A4","Plot G0 residuals",150,10,800,400);
TCanvas  *A5 = new TCanvas("A5","Plot G1 fit",150,10,800,400);
TCanvas  *A6 = new TCanvas("A6","Plot G1 residuals",150,10,800,400);

  
  

  //file name
  sprintf(filename,"%sBP_histos_M%s.root",folder_path.c_str(), module_str.c_str());
  cout << "Loading file " << filename << endl;
  TFile* f = new TFile((const char *)(filename),"READ");
   

  for (int j = 0; j < 245; j++) {
    TH2F* hist0=(TH2F*)f->Get(Form("avg_adcG0_map_%d",j));
    G0_pixel[j]=hist0->GetBinContent((column+1),(row+1));
    TH2F* hist0er=(TH2F*)f->Get(Form("avg_adcG0er_map_%d",j));
    G0_err_pixel[j]=hist0er->GetBinContent((column+1),(row+1));
    TH2F* hist1=(TH2F*)f->Get(Form("avg_adcG1_map_%d",j));
    G1_pixel[j]=hist1->GetBinContent((column+1),(row+1));
    TH2F* hist1er=(TH2F*)f->Get(Form("avg_adcG1er_map_%d",j));    
    G1_err_pixel[j]=hist1er->GetBinContent((column+1),(row+1));
    //cout << "Data for pixel "<< pixel << "is loaded" << endl;
    }  

  
  A2->SetLeftMargin(0.13);
  A2->SetRightMargin(0.05);

    vector<double> r0_adc;
    vector<double> r0_filter;

    vector<double> r0_adcerr;
    vector<double> r0_ferr;

    vector<double> r1_adc;
    vector<double> r1_filter;
    vector<double> r1_adcerr;
    vector<double> r1_ferr;

    for (int j = 0; j < 245; j++) {

      if (G0_pixel[j] != 0) {
	r0_filter.push_back(xs[j]);
	r0_ferr.push_back(0.);
	r0_adc.push_back(G0_pixel[j]);
	r0_adcerr.push_back(G0_err_pixel[j]);
      }	
      if (G1_pixel[j] != 0) {
	r1_filter.push_back(xs[j]);
	r1_ferr.push_back(0.);
	r1_adc.push_back(G1_pixel[j]);
	r1_adcerr.push_back(G1_err_pixel[j]);
      }
    cout << "Filter array is loaded" << endl;	
    }

    //TGraphErrors *grap_g0 = 0;
    //TGraphErrors *grap_g1 = 0;

    //TF1 *fit_g0 = 0;
    //TF1 *fit_g1 = 0;

    double rangemin_g0 = 0;
    double rangemax_g0 = 0;
    double rangemin_g1 = 0;
    double rangemax_g1 = 0;
     
    A2->cd(); 

    // define graphs
    if (r0_adc.size() > 1) {
      grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_ferr[0]),&(r0_adcerr[0]));
      grap_g0->SetMarkerStyle(20);
      grap_g0->SetMarkerColor(kBlue);
      grap_g0->SetLineColor(kBlue);

    }

    if (r1_adc.size() > 1) {
      grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_ferr[0]),&(r1_adcerr[0]));
      grap_g1->SetMarkerStyle(20);
      grap_g1->SetMarkerColor(kGreen+2);
      grap_g1->SetLineColor(kGreen+2);
    }

    //debug 
    //   cout << i <<" r0size= "<<  r0_adc.size()<< " r1size= "<<  r1_adc.size()<< endl;
    // plot the datapoints
    if (r1_adc.size() > 1) {

    grap_g1->GetXaxis()->SetTitle("Signal generator voltage [mV]");
    grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
    grap_g1->GetYaxis()->SetTitleOffset(0.9);
    grap_g1->SetMinimum(1000);
    grap_g1->SetMaximum(15000);
    grap_g1->GetXaxis()->SetLimits(0,7200);
    grap_g1->Draw("AP");
    if (r0_adc.size() > 1) {
           grap_g0->Draw("P");
           A2->Update();
    }
    
    //cout << "Number of point in g1 =" << r1_adc.size() << endl;
	//sprintf(savename,"plots/M%s/%s/pixel_%s_%d_M%s.png", module_str.c_str(),  plotfolder_str.c_str(),pixel_type.c_str(), i,module_str.c_str());
	//mapcanvas->SaveAs((const char *)(savename));

      }
    
    

    // define fit ranges and fit
    if (r0_adc.size() > 1) {
      A3->cd();
      
      rangemin_g0 = *min_element(r0_filter.begin(),r0_filter.end());
      rangemax_g0 = highestPointBeforeSwitching(r0_filter,r1_filter);

      if (rangemax_g0 > rangemin_g0) {
      
	grap_g0->Fit("pol1","QRC","",rangemin_g0,rangemax_g0);
	
    fit_g0 = (TF1*) grap_g0->GetFunction("pol1");
    
	if (fit_g0) {
      
	  fit_g0->SetLineColor(kBlue);
	  
      fit_g0->SetParName(0,"G0 const");
	  fit_g0->SetParName(1,"G0 grad");      
	  grap_g0->GetXaxis()->SetTitle("Signal generator voltage [mV]");
	  grap_g0->GetYaxis()->SetTitle("ADC [ADU]");
	  grap_g0->GetYaxis()->SetTitleOffset(0.9);
	  grap_g0->SetMinimum(1000);
	  grap_g0->SetMaximum(15000);
	  grap_g0->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));

              
	  // TPaveStats *st0 = (TPaveStats*)grap_g0->FindObject("stats");
	  // st0->SetX1NDC(0.2);
	  // st0->SetX2NDC(0.54);
	  // st0->SetY1NDC(0.18);
	  // st0->SetY2NDC(0.37);
	  // st0->SetBorderSize(0);
	  // st0->SetTextColor(kBlue);

      grap_g0->Draw("AP");
      fit_g0->Draw("same");
      A3->Update();
      
      //PlotCanvas();          
	  //sprintf(savename,"plots/M%s/%s/pixel_%s_%d_g0_M%s.png", module_str.c_str(),plotfolder_str.c_str(),  pixel_type.c_str(), i,module_str.c_str());
	  //mapcanvas->SaveAs((const char *)(savename));
      A4->cd();

	    vector<double> r0_adc_norm;
	    for (size_t j = 0; j < r0_adc.size(); j++) {
	      r0_adc_norm.push_back(r0_adc[j] - fit_g0->Eval(r0_filter[j]));
	    }

	    TGraphErrors *norm_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc_norm[0]),&(r0_ferr[0]),&(r0_adcerr[0]));
	    norm_g0->SetMarkerColor(kBlue);
	    norm_g0->SetLineColor(kBlue);

	    TF1* flat_g0 = new TF1("flat_g0","0",rangemin_g0,rangemax_g0);
	    flat_g0->SetLineColor(kBlue);

	    TF1* lin_g0_p1pc = new TF1("lin_g0_p1pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_p1pc->SetParameter(0,fit_g0->GetParameter(0)/100.);
	    lin_g0_p1pc->SetParameter(1,fit_g0->GetParameter(1)/100.);
	    lin_g0_p1pc->SetLineColor(kRed);

	    TF1* lin_g0_p05pc = new TF1("lin_g0_p05pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_p05pc->SetParameter(0,fit_g0->GetParameter(0)/200.);
	    lin_g0_p05pc->SetParameter(1,fit_g0->GetParameter(1)/200.);
	    lin_g0_p05pc->SetLineColor(kOrange+1);

	    TF1* lin_g0_p02pc = new TF1("lin_g0_p02pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_p02pc->SetParameter(0,fit_g0->GetParameter(0)/500.);
	    lin_g0_p02pc->SetParameter(1,fit_g0->GetParameter(1)/500.);
	    lin_g0_p02pc->SetLineColor(kOrange);

	    TF1* lin_g0_m1pc = new TF1("lin_g0_m1pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_m1pc->SetParameter(0,fit_g0->GetParameter(0)/-100.);
	    lin_g0_m1pc->SetParameter(1,fit_g0->GetParameter(1)/-100.);
	    lin_g0_m1pc->SetLineColor(kRed);

	    TF1* lin_g0_m05pc = new TF1("lin_g0_m05pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_m05pc->SetParameter(0,fit_g0->GetParameter(0)/-200.);
	    lin_g0_m05pc->SetParameter(1,fit_g0->GetParameter(1)/-200.);
	    lin_g0_m05pc->SetLineColor(kOrange+1);

	    TF1* lin_g0_m02pc = new TF1("lin_g0_m02pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
	    lin_g0_m02pc->SetParameter(0,fit_g0->GetParameter(0)/-500.);
	    lin_g0_m02pc->SetParameter(1,fit_g0->GetParameter(1)/-500.);
	    lin_g0_m02pc->SetLineColor(kOrange);

	    norm_g0->GetXaxis()->SetRangeUser(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
	  
	    norm_g0->SetMinimum(1.5*lin_g0_m02pc->Eval(rangemax_g0));
	    norm_g0->SetMaximum(1.5*lin_g0_p02pc->Eval(rangemax_g0));
	    
	    if (isHGX) {
	      norm_g0->SetMinimum(3*lin_g0_m02pc->Eval(rangemax_g0));
	      norm_g0->SetMaximum(3*lin_g0_p02pc->Eval(rangemax_g0));
	    }
	    
	    norm_g0->GetXaxis()->SetTitle("Signal generator voltage [mV]");
	    norm_g0->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	    norm_g0->GetYaxis()->SetTitleOffset(0.9);
        norm_g0->Draw("AP");
        flat_g0->Draw("same");
        lin_g0_p1pc->Draw("same");
        lin_g0_p05pc->Draw("same");
        lin_g0_p02pc->Draw("same");
        lin_g0_m1pc->Draw("same");
        lin_g0_m05pc->Draw("same");
        lin_g0_m02pc->Draw("same");
        norm_g0->Draw("P");
        A4->Update();

	    //sprintf(savename,"plots/M%s/%s/pixel_%s_%d_g0norm_M%s.png", module_str.c_str(), plotfolder_str.c_str(),pixel_type.c_str(), i,  module_str.c_str());
	    //mapcanvas->SaveAs((const char *)(savename));

	    //delete norm_g0;
	  
	 }     
    }
   }
   A5->cd();
    if (r1_adc.size() > 1) {

      rangemin_g1 = lowestPointAfterSwitching(r1_filter,r0_filter);
      rangemax_g1 = *max_element(r1_filter.begin(),r1_filter.end());
      rangemax_g1 = checkRangeMaxForAmplifierPlateau(rangemax_g1);

      if (rangemax_g1 > rangemin_g1) {

	grap_g1->Fit("pol1","QRC","",rangemin_g1,rangemax_g1);
	fit_g1 = (TF1*) grap_g1->GetFunction("pol1");
	if (fit_g1) {
	  fit_g1->SetLineColor(kGreen+2);
	  fit_g1->SetParName(0,"G1 const");
	  fit_g1->SetParName(1,"G1 grad");

      grap_g1->GetXaxis()->SetTitle("Signal generator voltage [mV]");
      grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
      grap_g1->GetYaxis()->SetTitleOffset(0.9);
      grap_g1->SetMinimum(1000);
      grap_g1->SetMaximum(15000);
      grap_g1->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
          

   //    TPaveStats *st0 = (TPaveStats*)grap_g1->FindObject("stats");
   //    st0->SetX1NDC(0.6);
   //    st0->SetX2NDC(0.94);
	  // st0->SetY1NDC(0.18);
	  // st0->SetY2NDC(0.37);
	  // st0->SetBorderSize(0);
	  // st0->SetTextColor(kGreen+2);
	  grap_g1->Draw("AP");
      fit_g1->Draw("same");
      A5->Update();
	  //sprintf(savename,"plots/M%s/%s/pixel_%s_%d_g1_M%s.png", module_str.c_str(), plotfolder_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  //mapcanvas->SaveAs((const char *)(savename));

	  vector<double> r1_adc_norm;
	  for (size_t j = 0; j < r1_adc.size(); j++) {
	    r1_adc_norm.push_back(r1_adc[j] - fit_g1->Eval(r1_filter[j]));
	  }
      A6->cd();
	  TGraphErrors *norm_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_ferr[0]),&(r1_adcerr[0]));
	  norm_g1->SetMarkerColor(kGreen+2);
	  norm_g1->SetLineColor(kGreen+2);

	  TF1* flat_g1 = new TF1("flat_gi","0",rangemin_g1,rangemax_g1);
	  flat_g1->SetLineColor(kGreen+2);

	  TF1* lin_g1_p02pc = new TF1("lin_g1_p02pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_p02pc->SetParameter(0,fit_g1->GetParameter(0)/500.);
	  lin_g1_p02pc->SetParameter(1,fit_g1->GetParameter(1)/500.);
	  lin_g1_p02pc->SetLineColor(kOrange);
      
          TF1* lin_g1_p01pc = new TF1("lin_g1_p01pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_p01pc->SetParameter(0,fit_g1->GetParameter(0)/1000.);
          lin_g1_p01pc->SetParameter(1,fit_g1->GetParameter(1)/1000.);
	  lin_g1_p01pc->SetLineColor(kYellow);

	  TF1* lin_g1_m02pc = new TF1("lin_g1_m02pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_m02pc->SetParameter(0,fit_g1->GetParameter(0)/-500.);
	  lin_g1_m02pc->SetParameter(1,fit_g1->GetParameter(1)/-500.);
	  lin_g1_m02pc->SetLineColor(kOrange);

	  TF1* lin_g1_m01pc = new TF1("lin_g1_m01pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_m01pc->SetParameter(0,fit_g1->GetParameter(0)/-1000.);
	  lin_g1_m01pc->SetParameter(1,fit_g1->GetParameter(1)/-1000.);
	  lin_g1_m01pc->SetLineColor(kYellow);

	  norm_g1->GetXaxis()->SetRangeUser(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
	  norm_g1->SetMinimum(1.5*lin_g1_m02pc->Eval(rangemin_g1));
	  norm_g1->SetMaximum(1.5*lin_g1_p02pc->Eval(rangemin_g1));
	  norm_g1->GetXaxis()->SetTitle("Signal generator voltage [mV]");
	  norm_g1->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	  norm_g1->GetYaxis()->SetTitleOffset(0.9);
      norm_g1->Draw("AP");
      flat_g1->Draw("same");
      lin_g1_p02pc->Draw("same");
      lin_g1_p01pc->Draw("same");
      lin_g1_m02pc->Draw("same");
      lin_g1_m01pc->Draw("same");

      norm_g1->Draw("P");


      A6->Update();
	  //sprintf(savename,"plots/M%s/%s/pixel_%s_%d_g1norm_M%s.png", module_str.c_str(),plotfolder_str.c_str(), pixel_type.c_str(), i,  module_str.c_str());
	  //mapcanvas->SaveAs((const char *)(savename));

        //delete norm_g1;  
	}
   }
}

    // get ratio measurements
    // if (fit_g0 && fit_g1) {

    //   double this_g0overg1 = fit_g0->GetParameter(1)/fit_g1->GetParameter(1);
    //   double this_g0overg1er = sqrt(pow(fit_g0->GetParError(1)/fit_g0->GetParameter(1),2) + pow(fit_g1->GetParError(1)/fit_g1->GetParameter(1),2));
    //   cout << "G0overG1 =" << this_g0overg1 << "+/-" << this_g0overg1er << endl;
    // } 
      
    rootapp.Run();
    return 0;

}