// calculate pedestal correction of fluo data for a given storage cell, only G0
// (can be called multiple times to handle multi-storage-cell data, each time with different Sc number as argument, i.e. supports manual multi-threading on the same data set)
// make correction and save spectrum per pixel
// for fitting call CuFluo_analysis_sc_singlethread_fits (with correct SC number) afterward

#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 <vector>

int main(int argc, char* argv[]) {

  jungfrauStyle();
  gStyle->SetOptFit(11);

  if (argc != 8) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify HG0 or G0" << endl;
    cout << "arg 3: specify data location folder" << endl;
    cout << "arg 4: specify pede file prefix" << endl;
    cout << "arg 5: specify data file prefix" << endl;
    cout << "arg 6: specify location of analysis root files" << endl;
    cout << "arg 7: specify storage cell number" << endl;
    cout << " " << endl;
    exit(1);
  }

  string module_str = argv[1];
  string gain_str = argv[2];
  bool isJF11=false;
  if (gain_str == "HG0JF11") {
    gain_str = "HG0";
    isJF11=true;
    }

  string data_loc = argv[3];
  string pede_file = argv[4];
  string data_file = argv[5];
  string anadata_loc = argv[6];

  bool save_to_afs = false;

  int createHistoFile = 1;

  char savename[256];
  char plotfolder[256];
  char rootdatafolder[256];
  int filen = 352;
  int pedefilen = 5;
  int this_storagecell = stoi(argv[7]); // needs compiler flag -std=c++11

  // 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/CuFluo
  sprintf(savename,"plots/M%s/CuFluo", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/CuFluo/G0
  sprintf(savename,"plots/M%s/CuFluo/G0", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/CuFluo/HG0
  sprintf(savename,"plots/M%s/CuFluo/HG0", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  sprintf(savename,"%s", anadata_loc.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

  if ( save_to_afs ) {

    cout << "Saving to afs." << endl;

    sprintf( plotfolder, "plots/M%s/CuFluo/%s", module_str.c_str(), gain_str.c_str() );
    sprintf( rootdatafolder, "data/M%s", module_str.c_str() );

  } else {

    cout << "Not saving to afs." << endl;

    // data/Mxxx
    sprintf( savename,"%s/data", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx
    sprintf( savename,"%s/plots", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo
    sprintf( savename,"%s/plots/CuFluo", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo/G0
    sprintf( savename,"%s/plots/CuFluo/G0", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo/HG0
    sprintf( savename,"%s/plots/CuFluo/HG0", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    sprintf( plotfolder, "%s/plots/CuFluo/%s", anadata_loc.c_str(), gain_str.c_str() );
    sprintf( rootdatafolder, "%s/data", anadata_loc.c_str() );

  }

  jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
  bool pixel_mask [NCH];
  pixelMaskObject->initialisePixelMask(pixel_mask);

  TCanvas *mapcanvas = new TCanvas("mapcanvas","",150,10,800,400);
  mapcanvas->SetLeftMargin(0.1);
  mapcanvas->SetRightMargin(0.13);
  mapcanvas->SetTopMargin(0.08);
  mapcanvas->SetBottomMargin(0.15);

  /*****************************************************************************************************************/
  /* NOTE: Since arrays, as they were used before, tend to cause memory leaks if they grow too large on the stack, */
  /*       it might be smarter to realize all array objects with std::vector instead.                              */
  /*****************************************************************************************************************/

  if (createHistoFile) {

    jungfrauFile *thisfile = new jungfrauFile();

    jungfrauPedestal* pedestalObject_SC = new jungfrauPedestal();
    pedestalObject_SC->pedestalSetNFrames(100);

    vector<uint16_t> pedestals16_G0_start( NCH ); //Track pedestal shifting over the course of data taking
    vector<double> pedeRMS16_G0( NCH );

    //cout << "Last element of pedeRMS16_G0 = " << pedeRMS16_G0[NCH-1] << endl;

    //double *pedeRMSary = pedeRMS16_G0.data();

    //cout << "Last element of pedeRMSary = " << pedeRMSary[NCH-1] << endl;

    { //scope to limit the following 2 histograms
      TH2F* pedestalsG0;
      TH2F* pedeRMSG0;
	
      Char_t *pedehistoname = new Char_t[50];
      snprintf( pedehistoname, 50, "pedestalsG0_sc%d", this_storagecell );
      pedestalsG0 = new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
      snprintf( pedehistoname, 50, "pedeRMSG0_sc%d", this_storagecell );
      pedeRMSG0 = new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

      for (int pedefilei = 0; pedefilei < pedefilen; ++pedefilei) {
    
	// open pede file
	sprintf( savename, "%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str() );
	thisfile->open( (char*)savename, pedefilei );
	
	//count events in file
	int nevents = 0;
	while ( thisfile->readNextFrame() ) {
	  nevents++;
	}
	thisfile->rewind();
	cout << "read " << nevents << " events" << endl;
	
	while ( thisfile->readNextFrame() ) {
	  // calculate pixel mask
	  pixelMaskObject->maskIfGainNot( 0, thisfile->getFrameDataHandle(), pixel_mask );
	  // caluclate pedestals
	  if( thisfile->currentSCnumber() == (uint64_t)this_storagecell ) {
	    pedestalObject_SC->addFrameToPedestalCalculation( thisfile->getFrameDataHandle() );
	  }
	}
	
	thisfile->close();
      
      } // end of loops over files
      
      for (int i = 0; i < NCH; ++i) {
	if (pixel_mask[i] == true) {
	  pedestalsG0->Fill( i%NC, i/NC, pedestalObject_SC->pedestalOfChannel(i) );
	  pedeRMSG0->Fill( i%NC, i/NC, pedestalObject_SC->rmsOfChannel(i) );
	  }
      }
      pedestalObject_SC->pedestalData((uint16_t*)(pedestals16_G0_start.data())); //this loads the pedestal data into the array pedestals16_G0_start
      pedestalObject_SC->pedestalRMSData(pedeRMS16_G0.data()); //same here
      pedestalObject_SC->pedestalResetUpdates();
      //pedestalObject[sci]->pedestalClear(); //I don't need to clear if I only calculate pedestal once.

      sprintf( savename, "%s/pixelmask_%s_M%s.png", plotfolder, gain_str.c_str(), module_str.c_str() );
      pixelMaskObject->plotPixelMask( pixel_mask, savename );
      cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
      mapcanvas->cd();
	
      pedestalsG0->GetXaxis()->SetTitle("Column");
      pedestalsG0->GetYaxis()->SetTitle("Row");
      pedestalsG0->GetYaxis()->SetTitleOffset(0.7);
      pedestalsG0->Draw("colz");
      sprintf(savename,"%s/pedeG0_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str());
      mapcanvas->SaveAs((const char *)(savename));

      pedestalsG0->Delete(); //not needed hereafter (try to save memory)
	
      pedeRMSG0->GetXaxis()->SetTitle("Column");
      pedeRMSG0->GetYaxis()->SetTitle("Row");
      pedeRMSG0->GetYaxis()->SetTitleOffset(0.7);
      pedeRMSG0->GetZaxis()->SetRangeUser(0,30);
      pedeRMSG0->Draw("colz");
      sprintf(savename,"%s/pedeRMSG0_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str());
      mapcanvas->SaveAs((const char *)(savename));

      pedeRMSG0->Delete(); //not needed hereafter (try to save memory)

      //Make sure histograms are deleted

    } // end of scope to limit pedestalsG0 and pedeRMSG0
    
    int adc2d_nbin= 1200;
    if (isJF11) adc2d_nbin= 1600;
    
    vector<TH2I*> adc2d; //declare
    adc2d.reserve(8);

    cout << "Initializing adc2d vector" << endl;
    
    //initialize
    for ( int adci = 0; adci < 8; ++adci ) {
      Char_t *_histoname = new Char_t[50];
      snprintf( _histoname, 50, "adc2d_%d_sc%d", adci + 1, this_storagecell );
      adc2d.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*adci-0.5 ), ( 65536*( adci+1 )-0.5 ) ) );
    }
    
    cout << "Done." << endl;
    
    //declare
    TH1D* adcpc_spec;
    TH2F* pede_updates;
    TH2F* pede_diff;
    
    //initialize
    Char_t *_histoname = new Char_t[50];
    snprintf( _histoname, 50, "adcpc_spec_sc%d", this_storagecell );
    adcpc_spec = new TH1D( _histoname, "", 300, 0, 3000 ); //spectrum for every sc, every 10000 frames
    snprintf( _histoname, 50, "pede_updates_sc%d", this_storagecell );
    pede_updates = new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
    snprintf( _histoname, 50, "pede_diff_sc%d", this_storagecell );
    pede_diff = new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

    for (int filei = 0; filei < filen; filei++) {

      // open data file
      sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str());
      thisfile->open((char*)savename, filei);
      
      while (thisfile->readNextFrame()) {
	
	uint16_t* imagedptr = thisfile->getFrameDataHandle();
	uint64_t scnumber = thisfile->currentSCnumber();
	
	if ( scnumber == (uint64_t)this_storagecell ) {
	  pedestalObject_SC->addG0FrameToPedestalCalculationWThreshold( imagedptr, pedestalObject_SC, pedeRMS16_G0.data() );

	  for (int i = 0; i < NCH; i++) {
	    
	    uint16_t gain = (imagedptr[i]&0xc000) >> 14;
	    
	    if (gain == 0) {
	      
	      int adcpc = (imagedptr[i]&0x3fff) - pedestalObject_SC->pedestalOfChannel(i);
	      adcpc_spec->Fill(adcpc);
	      
	      if (i < (65536*1)) {
		adc2d[0]->Fill(adcpc,i);
	      } else if (i < (65536*2)) {
		adc2d[1]->Fill(adcpc,i);
	      } else if (i < (65536*3)) {
		adc2d[2]->Fill(adcpc,i);
	      } else if (i < (65536*4)) {
		adc2d[3]->Fill(adcpc,i);
	      } else if (i < (65536*5)) {
		adc2d[4]->Fill(adcpc,i);
	      } else if (i < (65536*6)) {
		adc2d[5]->Fill(adcpc,i);
	      } else if (i < (65536*7)) {
		adc2d[6]->Fill(adcpc,i);
	      } else if (i < (65536*8)) {
		adc2d[7]->Fill(adcpc,i);
	      }
	    }
	  }
	}
      }
      thisfile->close();

      if ( filei%16 == 0 ) { // every 16th file equals every 10000 frames for each storage cell

	adcpc_spec->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
	adcpc_spec->Draw();
	mapcanvas->SetLogy();
	sprintf(savename,"%s/adcpc_spec_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	mapcanvas->SaveAs((const char *)(savename));
	mapcanvas->SetLogy(0);
	adcpc_spec->Reset();

	pede_updates->Reset();
	pede_diff->Reset();
	for (int i = 0; i < NCH; i++) {
	  pede_updates->Fill( i%NC, i/NC, pedestalObject_SC->pedestalUpdates(i) );
	  pede_diff->Fill( i%NC, i/NC, pedestalObject_SC->pedestalOfChannel(i) - pedestals16_G0_start[i] );
	  pedestals16_G0_start[i] = pedestalObject_SC->pedestalOfChannel(i);
	}

	pede_updates->GetXaxis()->SetTitle("Column");
	pede_updates->GetYaxis()->SetTitle("Row");
	pede_updates->GetYaxis()->SetTitleOffset(0.7);
	pede_updates->GetZaxis()->SetRangeUser(0,10000);
	pede_updates->Draw("colz");
	sprintf(savename,"%s/pede_updates_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	mapcanvas->SaveAs((const char *)(savename));

	pede_diff->GetXaxis()->SetTitle("Column");
	pede_diff->GetYaxis()->SetTitle("Row");
	pede_diff->GetYaxis()->SetTitleOffset(0.7);
	pede_diff->GetZaxis()->SetRangeUser(-40,40);
	pede_diff->Draw("colz");
	sprintf(savename,"%s/pede_diff_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	mapcanvas->SaveAs((const char *)(savename));
	  
	pedestalObject_SC->pedestalResetUpdates();
	
      
      } //end of if filei%16

    } // end of file loop

    sprintf(savename,"%s/CuFluo_%s_sc%d_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), this_storagecell, filen-1);
    TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
    for ( int adci = 0; adci < 8; ++adci ) {
	adc2d[adci]->Write();

	//In theory, I believe, adc2d can now also be destroyed since it has already been written to file
	adc2d[adci]->Delete();
    }
    //All pede updates have been written, so they can be destroyed.
    pede_updates->Delete();
    pede_diff->Delete();
      
    saved_file->Close();
    
  } // end if creatHistoFile


}