// 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/jungfrauPixelMask.C"
#include "sls_detector_calibration/jungfrauPedestal.C"

#include "sls_detector_calibration/energyCalibration.h"
#include "sls_detector_calibration/energyCalibration.cpp"

#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"
#include "TLegend.h"
#include "TPaveText.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;


//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 != 5) {  
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: gain settings" << endl;
    cout << "arg 3: specify pixel x position" << endl;
    cout << "arg 4: specify pixel y position" << endl;
    exit(1);
  } //uncomment for VH 210906
  string module_str = argv[1];
  string gain_str = argv[2];

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


  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; 

  bool isJF11=false;
  if (gain_str == "HG0JF11") {
    gain_str = "HG0";
    isJF11=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"; 
  }
 
 
  //char savename[128]; //uncomment for SR
  //char filename[128]; //uncomment for SR
  char filename[256]; //uncomment for VH 210902


  int low_ADU_peak = 0;
  int high_ADU_peak = 0;
  if (gain_str == "HG0") {
    
    
    low_ADU_peak = 700;
    high_ADU_peak = 900;
    if  (isJF11) {
        low_ADU_peak = 850;
    high_ADU_peak = 1350;

    }

  } else if (gain_str == "G0") {
    low_ADU_peak = 250;
    high_ADU_peak = 400;
  }

  int low_bin_noise = 101;
  int high_bin_noise = 301;
  int low_bin_peak = 0;
  int high_bin_peak = 0;
  if (gain_str == "HG0") {
    low_bin_peak = 701;
    high_bin_peak = 1200;
     if  (isJF11) {
        low_bin_peak = 801;
	high_bin_peak = 1451;

    }

} else if (gain_str == "G0") {
    low_bin_peak = 301;
    high_bin_peak = 651;
  }

TCanvas  *A2 = new TCanvas("A2","Plot histogram",150,10,800,400);
TCanvas  *A3 = new TCanvas("A3","Plot noise",150,10,800,400);

 A2->SetLeftMargin(0.1);
 A2->SetRightMargin(0.13);
 A2->SetTopMargin(0.08);
 A2->SetBottomMargin(0.15);

 A3->SetLeftMargin(0.1);
 A3->SetRightMargin(0.13);
 A3->SetTopMargin(0.08);
 A3->SetBottomMargin(0.15);

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

 int chip;
 int CH = row*1024+column;
  
 if (CH < (65536*1)) {
     chip = 1;     
 } else if (CH < (65536*2)) {
     chip = 2;
 } else if (CH < (65536*3)) {
     chip = 3;    
 } else if (CH < (65536*4)) {     
     chip = 4;  
 } else if (CH < (65536*5)) {
     chip = 5;
 } else if (CH < (65536*6)) {
     chip = 6;  
 } else if (CH < (65536*7)) {
     chip = 7;   
 } else if (CH < (65536*8)){
     chip = 8;
 }
  cout<< "Chip "<< chip << " Channel number " << CH << endl;  

 TH2I* hist0=(TH2I*)f->Get(Form("adc2d_%d",chip));
 cout << "Creating histogram for " << Form("adc2d_%d",chip) << endl;
 TH1D* proj = hist0->ProjectionX("bin1",CH-(65536*(chip-1))+1,CH-(65536*(chip-1))+1);
 cout << "Data for pixel "<< pixel << " is loaded" << endl;


if (proj->Integral(low_bin_noise,high_bin_noise)!=0 && proj->Integral(low_bin_peak,high_bin_peak)!=0) {
	
    A2->cd();
	cout << "Canvas Noise" << endl;
    // noise
	TH1D *proj_noise = dynamic_cast<TH1D*>(proj->Rebin(4,"proj_noise"));
	proj_noise->SetStats(kTRUE);
	proj_noise->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_noise),proj->GetBinLowEdge(high_bin_noise+1));
	proj_noise->Fit("gaus","Q");
	TF1 *fit = proj_noise->GetFunction("gaus");

    proj_noise->Draw();	
	A2->Update();
    proj_noise->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
	proj_noise->GetXaxis()->SetRangeUser(-100,150);
	fit->SetParNames("N_{#gamma}", "Peak pos", "Noise RMS");
	TPaveStats *st0 = (TPaveStats*)proj_noise->FindObject("stats");
	st0->SetX1NDC(0.53);
	st0->SetX2NDC(0.94);
	st0->SetY1NDC(0.75);
	st0->SetY2NDC(0.94);
	st0->SetBorderSize(0);
	st0->SetTextSize(0.04);
	A2->Modified();
	A2->Update();

    A3->cd();

	// peak
	TH1D *proj_peak = dynamic_cast<TH1D*>(proj->Rebin(4,"proj_peak"));
	proj_peak->SetStats(kTRUE);
	proj_peak->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));


    Double_t mypar[8];
	mypar[0] = 0.0;
	mypar[1] = 0.0;
	mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
	if (gain_str == "G0") {
	  mypar[3] = 16.;
	} else if (gain_str == "HG0") {
	  mypar[3] = 29.;
	}
	mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
	if (gain_str == "G0") {
	  mypar[5] = 0.17;
	} else if (gain_str == "HG0") {
	  mypar[5] = 0.14;
	}
	mypar[6] = 1.12;
	if (gain_str == "G0") {
	  mypar[7] = 0.12;
	} else if (gain_str == "HG0") {
	  mypar[7] = 0.14;
	}

    Double_t emypar[8];
	energyCalibration *thiscalibration = new energyCalibration();
	thiscalibration->setScanSign(1);
	thiscalibration->setStartParametersKb(mypar);
	thiscalibration->fixParameter(0,0.); // no background
	thiscalibration->fixParameter(1,0.);
	TF1* fittedfun = thiscalibration->fitSpectrumKb(proj_peak,mypar,emypar);

	TF1 *gaus_Ka = new TF1("gaus_Ka","gaus",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
	gaus_Ka->SetParameters(mypar[4],mypar[2],mypar[3]);
	gaus_Ka->SetLineColor(kBlue);
	
	TF1 *erfc_Ka = new TF1("erfc_Ka","[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
	erfc_Ka->SetParameters(mypar[4]*mypar[5], 1, mypar[2], mypar[3]);
	erfc_Ka->SetLineColor(kOrange);

	TF1 *gaus_Kb = new TF1("gaus_Kb","gaus",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
	gaus_Kb->SetParameters(mypar[4]*mypar[7],mypar[6]*mypar[2],mypar[3]);
	gaus_Kb->SetLineColor(kGreen+2);

	TF1 *erfc_Kb = new TF1("erfc_Kb","[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
	erfc_Kb->SetParameters(mypar[4]*mypar[7]*mypar[5], 1, mypar[6]*mypar[2], mypar[3]);
	erfc_Kb->SetLineColor(kOrange+7);

	proj_peak->Draw();
	A3->Update();
	erfc_Kb->Draw("same");
	erfc_Ka->Draw("same");
	gaus_Kb->Draw("same");
	gaus_Ka->Draw("same");
	fittedfun->Draw("same");
	A3->Update();
	proj_peak->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
	fittedfun->SetParNames("Bkg height", "Bkg grad", "K_{#alpha} pos", "Noise RMS", "K_{#alpha} height", "CS", "K_{#beta}/K_{#alpha} pos", "K_{#beta} frac");
	TPaveStats *st = (TPaveStats*)proj_peak->FindObject("stats");
	st->SetX1NDC(0.15);
	st->SetX2NDC(0.55);
	st->SetY1NDC(0.7);
	st->SetY2NDC(0.94);
	st->SetBorderSize(0);
	st->SetTextSize(0.04);
	A3->Modified();
	A3->Update();
} else {
	std::cout << "Masked pixel";
      }

    
rootapp.Run();
return 0;

}