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

#include "TGraphErrors.h"
#include "TF1.h"
#include "TPaveStats.h"

#include <sstream>

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

  jungfrauStyle();
  gROOT->SetBatch(1);
  gStyle->SetOptFit(11);

  if (argc != 5) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify data location" << endl;
    cout << "arg 3: specify scan id" << endl;
    cout << "arg 4: specify target pixel" << endl;
    cout << " " << endl;
    exit(1);
  }

  string module_str = argv[1];
  string data_file = argv[2];
  string scan_id = argv[3];

  int target_pixel = 0;
  stringstream target_pixel_ss(argv[4]);
  target_pixel_ss >> target_pixel;

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

  bool saturation_flag = true;

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

  jungfrauPedestal *pedestalObjectG0 = new jungfrauPedestal();
  pedestalObjectG0->pedestalSetNFrames(25);
  jungfrauPedestal *pedestalObjectG1 = new jungfrauPedestal();
  pedestalObjectG1->pedestalSetNFrames(25);
  jungfrauPedestal *pedestalObjectG2 = new jungfrauPedestal();
  pedestalObjectG2->pedestalSetNFrames(25);

  jungfrauFile *thisfile = new jungfrauFile();
  
  TH2F* pedestalsG0 = new TH2F("pedestalsG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* pedestalsG1 = new TH2F("pedestalsG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* pedestalsG2 = new TH2F("pedestalsG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  TH2F* pedeRMSG0 = new TH2F("pedeRMSG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* pedeRMSG1 = new TH2F("pedeRMSG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* pedeRMSG2 = new TH2F("pedeRMSG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  TCanvas *c1 = new TCanvas("c1","");
  
  TCanvas *viscanvas = new TCanvas("viscanvas","");
  viscanvas->SetLeftMargin(0.15);
  viscanvas->SetRightMargin(0.2);
  viscanvas->SetTopMargin(0.08);
  viscanvas->SetBottomMargin(0.15);

  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);

  char savename[128];

  TH2F *adc_map = new TH2F("adc_map","",20,(target_pixel%1024)-10-0.5,(target_pixel%1024)+10-0.5,20,(target_pixel/1024)-10-0.5,(target_pixel/1024)+10-0.5);
  TH2F *gain_map = new TH2F("gain_map","",20,(target_pixel%1024)-10-0.5,(target_pixel%1024)+10-0.5,20,(target_pixel/1024)-10-0.5,(target_pixel/1024)+10-0.5);

  TH1F* g0overg1_hist = new TH1F("g0overg1_hist","",100,-34,-24);
  TH1F* g0overg1er_hist = new TH1F("g0overg1er_hist","",100,0,0.1);
  TH1F* g1overg2_hist = new TH1F("g1overg2_hist","",100,10,25);
  TH1F* g1overg2er_hist = new TH1F("g1overg2er_hist","",100,0,0.01);
    
  sprintf(savename,"%s_%%6.6d.dat", data_file.c_str());
  thisfile->open((char*)savename, 0);

  // total frame number check
  int expected_totalFrames = 0;
  if (module_str == "002" && (scan_id == "1" || scan_id == "2")) {
    expected_totalFrames = 1499;
  } else if (module_str == "002") {
    expected_totalFrames = 7499;
  } else {
    expected_totalFrames = 7500;
  }

  int totalFrames = 0;

  while(thisfile->readNextFrame()) {
    totalFrames++;
  }

  thisfile->rewind();

  if (totalFrames == expected_totalFrames) {
    cout << "expected number of frames " << totalFrames << endl;
  } else {
    cout << "unexpected number of frames " << totalFrames << endl;
    cout << "code expected frames " << expected_totalFrames << endl;
    cout << "exiting to make sure you notice" << endl;
    exit(1);
  }

  int framecounter = 0;

  Double_t adc_avg_g0[60];
  Double_t adc_avg_g1[60];
  Double_t adc_avg_g2[60];
  Double_t adcer_avg_g0[60];
  Double_t adcer_avg_g1[60];
  Double_t adcer_avg_g2[60];

  if (not (module_str == "002" && (scan_id == "1" || scan_id == "2"))) {

    for (int i = 0; i < 2000; i++) {
      thisfile->readNextFrame();
      pixelMaskObject->maskIfGainNot(0, thisfile->getFrameDataHandle(), pixel_mask);
    }
    cout << "after G0, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
    for (int i = 0; i < 2000; i++) {
      thisfile->readNextFrame();
      pixelMaskObject->maskIfGainNot(1, thisfile->getFrameDataHandle(), pixel_mask);
    }
    cout << "after G1, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
    for (int i = 0; i < 2000; i++) {
      thisfile->readNextFrame();
      pixelMaskObject->maskIfGainNot(3, thisfile->getFrameDataHandle(), pixel_mask);
    }
    cout << "after G2, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
    thisfile->rewind();
      
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pixelmask_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    pixelMaskObject->plotPixelMask(pixel_mask,savename);
    cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
    // calculate pedestals
    for (int j = 0; j < 2000; j++) {
      thisfile->readNextFrame();
      pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
    }
    for (int i = 0; i < NCH; i++) {
      if (pixel_mask[i] == true) {
	pedestalsG0->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
	pedeRMSG0->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
      }
    }
    pedestalObject->pedestalClear();
    
    for (int j = 0; j < 2000; j++) {
      thisfile->readNextFrame();
      pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
    }
    for (int i = 0; i < NCH; i++) {
      if (pixel_mask[i] == true) {
	pedestalsG1->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
	pedeRMSG1->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
      }
    }
    pedestalObject->pedestalClear();

    for (int j = 0; j < 2000; j++) {
      thisfile->readNextFrame();
      pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
    }
    for (int i = 0; i < NCH; i++) {
      if (pixel_mask[i] == true) {
	pedestalsG2->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
	pedeRMSG2->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
      }
    }
    pedestalObject->pedestalClear();

    mapcanvas->cd();

    pedestalsG0->GetXaxis()->SetTitle("Column");
    pedestalsG0->GetYaxis()->SetTitle("Row");
    pedestalsG0->GetYaxis()->SetTitleOffset(0.7);
    pedestalsG0->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeG0_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));

    pedestalsG1->GetXaxis()->SetTitle("Column");
    pedestalsG1->GetYaxis()->SetTitle("Row");
    pedestalsG1->GetYaxis()->SetTitleOffset(0.7);
    pedestalsG1->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeG1_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));

    pedestalsG2->GetXaxis()->SetTitle("Column");
    pedestalsG2->GetYaxis()->SetTitle("Row");
    pedestalsG2->GetYaxis()->SetTitleOffset(0.7);
    pedestalsG2->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeG2_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));

    pedeRMSG0->GetXaxis()->SetTitle("Column");
    pedeRMSG0->GetYaxis()->SetTitle("Row");
    pedeRMSG0->GetYaxis()->SetTitleOffset(0.7);
    pedeRMSG0->GetZaxis()->SetRangeUser(0,25);
    pedeRMSG0->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeRMSG0_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));

    pedeRMSG1->GetXaxis()->SetTitle("Column");
    pedeRMSG1->GetYaxis()->SetTitle("Row");
    pedeRMSG1->GetYaxis()->SetTitleOffset(0.7);
    pedeRMSG1->GetZaxis()->SetRangeUser(0,15);
    pedeRMSG1->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeRMSG1_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));

    pedeRMSG2->GetXaxis()->SetTitle("Column");
    pedeRMSG2->GetYaxis()->SetTitle("Row");
    pedeRMSG2->GetYaxis()->SetTitleOffset(0.7);
    pedeRMSG2->GetZaxis()->SetRangeUser(0,15);
    pedeRMSG2->Draw("colz");
    sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/pedeRMSG2_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
    mapcanvas->SaveAs((const char *)(savename));
  }

  while (thisfile->readNextFrame()) {

    uint16_t* imagedptr = thisfile->getFrameDataHandle();
    uint16_t adc = imagedptr[target_pixel]&0x3fff;
    uint16_t gain = (imagedptr[target_pixel]&0xc000) >> 14;

    for (int i = 0; i < NCH; i++) {
      adc_map->Fill(i%1024,i/1024,imagedptr[i]&0x3fff);
      gain_map->Fill(i%1024,i/1024,((imagedptr[i]&0xc000) >> 14));
    }

    pedestalObjectG0->addG0FrameToPedestalCalculation(imagedptr);
    pedestalObjectG1->addG1FrameToPedestalCalculation(imagedptr);
    pedestalObjectG2->addG2FrameToPedestalCalculation(imagedptr);

    if (gain == 3 && adc == 0) {
      saturation_flag = false;
    }

    framecounter++;

    if (framecounter%25 == 0 || framecounter == expected_totalFrames) {

      // catch the last frame (ctrl-c exit)
      if (framecounter == 1499) {
	framecounter = 1500;
      } else if (framecounter == 7499) {
	framecounter = 7500;
      }

      viscanvas->cd();

      adc_map->Scale(1./25);
      gain_map->Scale(1./25);

      adc_map->GetXaxis()->SetTitle("Column");
      adc_map->GetYaxis()->SetTitle("Row");
      adc_map->Draw("colz");
      adc_map->GetZaxis()->SetRangeUser(0,16400);
      sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/perStep/map_adc_%i.png", module_str.c_str(), scan_id.c_str(), framecounter/25);
      viscanvas->SaveAs((const char *)(savename));

      gain_map->GetXaxis()->SetTitle("Column");
      gain_map->GetYaxis()->SetTitle("Row");
      gain_map->Draw("colz");
      gain_map->GetZaxis()->SetRangeUser(0,3);
      sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/perStep/gain_map_%i.png", module_str.c_str(), scan_id.c_str(), framecounter/25);
      viscanvas->SaveAs((const char *)(savename));

      adc_avg_g0[(framecounter/25)-1] = pedestalObjectG0->pedestalOfChannel(target_pixel);
      adcer_avg_g0[(framecounter/25)-1] = pedestalObjectG0->semOfChannel(target_pixel);
      adc_avg_g1[(framecounter/25)-1] = pedestalObjectG1->pedestalOfChannel(target_pixel);
      adcer_avg_g1[(framecounter/25)-1] = pedestalObjectG1->semOfChannel(target_pixel);
      if (saturation_flag == true) {
	adc_avg_g2[(framecounter/25)-1] = pedestalObjectG2->pedestalOfChannel(target_pixel);
	adcer_avg_g2[(framecounter/25)-1] = pedestalObjectG2->semOfChannel(target_pixel);
      }

      adc_map->Reset();
      gain_map->Reset();
 
      pedestalObjectG0->pedestalClear();
      pedestalObjectG1->pedestalClear();
      pedestalObjectG2->pedestalClear();

      saturation_flag = true;
    }
  }
  thisfile->close();

  // FilterCalibration_IR_April2014.txt
  // this one looks better than the later one
  Double_t filter[60];
  filter[0]=29;
  filter[1]=21;
  filter[2]=31;
  filter[3]=43;
  filter[4]=60;
  filter[5]=91;
  filter[6]=69;
  filter[7]=102;
  filter[8]=136;
  filter[9]=196;
  filter[10]=425;
  filter[11]=311;
  filter[12]=462;
  filter[13]=653;
  filter[14]=926;
  filter[15]=1423;
  filter[16]=1072;
  filter[17]=1592;
  filter[18]=2142;
  filter[19]=3085;
  filter[20]=729;
  filter[21]=533;
  filter[22]=793;
  filter[23]=1121;
  filter[24]=1588;
  filter[25]=2439;
  filter[26]=1842;
  filter[27]=2730;
  filter[28]=3663;
  filter[29]=5271;
  filter[30]=8102;
  filter[31]=5933;
  filter[32]=8789;
  filter[33]=12350;
  filter[34]=17358;
  filter[35]=26300;
  filter[36]=20029;
  filter[37]=29414;
  filter[38]=39202;
  filter[39]=55724;
  filter[40]=15697;
  filter[41]=11541;
  filter[42]=16976;
  filter[43]=23866;
  filter[44]=33478;
  filter[45]=50567;
  filter[46]=38552;
  filter[47]=56394;
  filter[48]=74897;
  filter[49]=106023;
  filter[50]=157384;
  filter[51]=117677;
  filter[52]=171101;
  filter[53]=236386;
  filter[54]=327248;
  filter[55]=492781;
  filter[56]=379641;
  filter[57]=546927;
  filter[58]=717203;
  filter[59]=1000000;

  // Calibration_IR_15April2016
  /*
    Double_t filter[60];
    filter[0]=42.578;
    filter[1]=31.397;
    filter[2]=46.218;
    filter[3]=64.466;
    filter[4]=90.267;
    filter[5]=136.928;
    filter[6]=104.231;
    filter[7]=152.697;
    filter[8]=202.840;
    filter[9]=290.595;
    filter[10]=616.125;
    filter[11]=454.778;
    filter[12]=669.518;
    filter[13]=933.511;
    filter[14]=1289.725;
    filter[15]=1945.277;
    filter[16]=1498.417;
    filter[17]=2204.015;
    filter[18]=2938.539;
    filter[19]=4208.522;
    filter[20]=1050.173;
    filter[21]=775.664;
    filter[22]=1142.794;
    filter[23]=1595.188;
    filter[24]=2231.033;
    filter[25]=3380.634;
    filter[26]=2575.221;
    filter[27]=3768.117;
    filter[28]=5003.264;
    filter[29]=7164.062;
    filter[30]=10875.166;
    filter[31]=8046.250;
    filter[32]=11845.280;
    filter[33]=16525.147;
    filter[34]=23061.385;
    filter[35]=34755.307;
    filter[36]=26580.837;
    filter[37]=38843.108;
    filter[38]=50913.202;
    filter[39]=74590.827;
    filter[40]=20862.210;
    filter[41]=15386.015;
    filter[42]=22403.187;
    filter[43]=32168.090;
    filter[44]=44968.315;
    filter[45]=68190.429;
    filter[46]=52051.379;
    filter[47]=76292.639;
    filter[48]=101288.429;
    filter[49]=145090.792;
    filter[50]=217266.543;
    filter[51]=160650.860;
    filter[52]=236126.004;
    filter[53]=329205.971;
    filter[54]=460252.653;
    filter[55]=696729.217;
    filter[56]=531860.655;
    filter[57]=777668.189;
    filter[58]=1030434.153;
    filter[59]=1471487.831;
  */

  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;

  vector<double> r2_adc;
  vector<double> r2_filter;
  vector<double> r2_adcerr;
  vector<double> r2_ferr;

  // filter out points at zero
  for (int i = 0; i < 60; i++) {

    double this_g0avg = adc_avg_g0[i];
    double this_g1avg = adc_avg_g1[i];
    double this_g2avg = adc_avg_g2[i];
    double this_g0avger = adcer_avg_g0[i];
    double this_g1avger = adcer_avg_g1[i];
    double this_g2avger = adcer_avg_g2[i];

    if (this_g0avg != 0) {
      r0_filter.push_back(filter[i]);
      r0_ferr.push_back(0.);
      r0_adc.push_back(this_g0avg);
      r0_adcerr.push_back(this_g0avger);
    }
    if (this_g1avg != 0) {
      r1_filter.push_back(filter[i]);
      r1_ferr.push_back(0.);
      r1_adc.push_back(this_g1avg);
      r1_adcerr.push_back(this_g1avger);
    }
    if (this_g2avg != 0) {
      r2_filter.push_back(filter[i]);
      r2_ferr.push_back(0.);
      r2_adc.push_back(this_g2avg);
      r2_adcerr.push_back(this_g2avger);
    }
  }

  // define graphs
  TGraphErrors *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);

  TGraphErrors *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);

  TGraphErrors *grap_g2 = new TGraphErrors(r2_adc.size(),&(r2_filter[0]),&(r2_adc[0]),&(r2_ferr[0]),&(r2_adcerr[0]));
  grap_g2->SetMarkerStyle(20);
  grap_g2->SetMarkerColor(kRed);
  grap_g2->SetLineColor(kRed);

  mapcanvas->cd();
  mapcanvas->SetLeftMargin(0.13);
  mapcanvas->SetRightMargin(0.05);

  mapcanvas->SetLogx();
  grap_g0->SetMinimum(0);
  grap_g0->SetMaximum(16000);
  grap_g0->GetXaxis()->SetTitle("Intensity [AU]");
  grap_g0->GetYaxis()->SetTitle("ADC [ADC]");
  grap_g0->GetYaxis()->SetTitleOffset(0.9);
  grap_g0->Draw("AP");
  grap_g1->Draw("P");
  grap_g2->Draw("P");
  grap_g0->GetXaxis()->SetLimits(5,2000000);
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));
  //mapcanvas->SetLogx(0);

  double rangemin_g0 = *min_element(r0_filter.begin(),r0_filter.end());
  double rangemax_g0 = highestPointBeforeSwitching(r0_filter,r1_filter);
  double rangemin_g1 = lowestPointAfterSwitching(r1_filter,r0_filter);
  double rangemax_g1 = highestPointBeforeSwitching(r1_filter,r2_filter);
  double rangemin_g2 = lowestPointAfterSwitching(r2_filter,r1_filter);
  double rangemax_g2 = *max_element(r2_filter.begin(),r2_filter.end());

  // fit
  grap_g0->Fit("pol1","QRC","",rangemin_g0,rangemax_g0);
  TF1 *fit_g0 = (TF1*) grap_g0->GetFunction("pol1");
  fit_g0->SetLineColor(kBlue);
  fit_g0->SetParName(0,"G0 const");
  fit_g0->SetParName(1,"G0 grad");

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

  grap_g2->Fit("pol1","QRC","",rangemin_g2,rangemax_g2);
  TF1 *fit_g2 = (TF1*) grap_g2->GetFunction("pol1");
  fit_g2->SetLineColor(kRed);
  fit_g2->SetParName(0,"G2 const");
  fit_g2->SetParName(1,"G2 grad");

  // draw
  grap_g0->GetXaxis()->SetTitle("Intensity [AU]");
  grap_g0->GetYaxis()->SetTitle("ADC [ADC]");
  grap_g0->GetYaxis()->SetTitleOffset(0.9);
  grap_g0->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
  grap_g0->Draw("AP");
  mapcanvas->Update();
  TPaveStats *st0 = (TPaveStats*)grap_g0->FindObject("stats");
  st0->SetX1NDC(0.6);
  st0->SetX2NDC(0.94);
  st0->SetY1NDC(0.17);
  st0->SetY2NDC(0.36);
  st0->SetBorderSize(0);
  st0->SetTextColor(kBlue);
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g0_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  grap_g1->GetXaxis()->SetTitle("Intensity [AU]");
  grap_g1->GetYaxis()->SetTitle("ADC [ADC]");
  grap_g1->GetYaxis()->SetTitleOffset(0.9);
  grap_g1->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
  grap_g1->Draw("AP");
  mapcanvas->Update();
  TPaveStats *st1 = (TPaveStats*)grap_g1->FindObject("stats");
  st1->SetX1NDC(0.2);
  st1->SetX2NDC(0.54);
  st1->SetY1NDC(0.17);
  st1->SetY2NDC(0.36);
  st1->SetBorderSize(0);
  st1->SetTextColor(kGreen+2);
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g1_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  grap_g2->GetXaxis()->SetTitle("Intensity [AU]");
  grap_g2->GetYaxis()->SetTitle("ADC [ADC]");
  grap_g2->GetYaxis()->SetTitleOffset(0.9);
  grap_g2->GetXaxis()->SetLimits(*min_element(r2_filter.begin(),r2_filter.end()),*max_element(r2_filter.begin(),r2_filter.end()));
  grap_g2->Draw("AP");
  mapcanvas->Update();
  TPaveStats *st2 = (TPaveStats*)grap_g2->FindObject("stats");
  st2->SetX1NDC(0.2);
  st2->SetX2NDC(0.54);
  st2->SetY1NDC(0.17);
  st2->SetY2NDC(0.36);
  st2->SetBorderSize(0);
  st2->SetTextColor(kRed);
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g2_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  // residuals
  vector<double> r0_adc_norm;
  vector<double> r1_adc_norm;
  vector<double> r2_adc_norm;

  for (size_t i = 0; i < r0_adc.size(); i++) {
    r0_adc_norm.push_back(r0_adc[i] - fit_g0->Eval(r0_filter[i]));
  }
  for (size_t i = 0; i < r1_adc.size(); i++) {
    r1_adc_norm.push_back(r1_adc[i] - fit_g1->Eval(r1_filter[i]));
  }
  for (size_t i = 0; i < r2_adc.size(); i++) {
    r2_adc_norm.push_back(r2_adc[i] - fit_g2->Eval(r2_filter[i]));
  }

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

  TGraphErrors *norm_g1 = new TGraphErrors(r1_adc_norm.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_ferr[0]),&(r1_adcerr[0]));
  norm_g1->SetMarkerStyle(20);
  norm_g1->SetMarkerColor(kGreen+2);
  norm_g1->SetLineColor(kGreen+2);

  TGraphErrors *norm_g2 = new TGraphErrors(r2_adc_norm.size(),&(r2_filter[0]),&(r2_adc_norm[0]),&(r2_ferr[0]),&(r2_adcerr[0]));
  norm_g2->SetMarkerStyle(20);
  norm_g2->SetMarkerColor(kRed);
  norm_g2->SetLineColor(kRed);

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

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

  TF1* flat_g2 = new TF1("flat_g2","0",rangemin_g2,rangemax_g2);
  flat_g2->SetLineColor(kRed);

  norm_g0->GetXaxis()->SetTitle("Intensity [AU]");
  norm_g0->GetYaxis()->SetTitle("Residual [ADC]");
  norm_g0->GetYaxis()->SetTitleOffset(0.9);
  norm_g0->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
  norm_g0->Draw("AP");
  flat_g0->Draw("same");
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g0_norm_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  norm_g1->GetXaxis()->SetTitle("Intensity [AU]");
  norm_g1->GetYaxis()->SetTitle("Residual [ADC]");
  norm_g1->GetYaxis()->SetTitleOffset(0.9);
  norm_g1->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
  norm_g1->Draw("AP");
  flat_g1->Draw("same");
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g1_norm_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  norm_g2->GetXaxis()->SetTitle("Intensity [AU]");
  norm_g2->GetYaxis()->SetTitle("Residual [ADC]");
  norm_g2->GetYaxis()->SetTitleOffset(0.9);
  norm_g2->GetXaxis()->SetLimits(*min_element(r2_filter.begin(),r2_filter.end()),*max_element(r2_filter.begin(),r2_filter.end()));
  norm_g2->Draw("AP");
  flat_g2->Draw("same");
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/intensity_adc_g2_norm_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  double this_g0overg1 = fit_g0->GetParameter(1)/fit_g1->GetParameter(1);
  double this_g0overg1er = sqrt((fit_g0->GetParError(1)*fit_g0->GetParError(1)) + (fit_g1->GetParError(1)*fit_g1->GetParError(1)));
  double this_g1overg2 = fit_g1->GetParameter(1)/fit_g2->GetParameter(1);
  double this_g1overg2er = sqrt((fit_g1->GetParError(1)*fit_g1->GetParError(1)) + (fit_g2->GetParError(1)*fit_g2->GetParError(1)));

  c1->cd();

  g0overg1_hist->Fill(this_g0overg1);
  g0overg1_hist->GetXaxis()->SetTitle("G0 / G1");
  g0overg1_hist->Draw();
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/g0overg1_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g0overg1er_hist->Fill(this_g0overg1er);
  g0overg1er_hist->GetXaxis()->SetTitle("#sigma_{G0/G1}");
  g0overg1er_hist->Draw();
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/g0overg1er_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g1overg2_hist->Fill(this_g1overg2);
  g1overg2_hist->GetXaxis()->SetTitle("G1 / G2");
  g1overg2_hist->Draw();
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/g1overg2_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));
    
  g1overg2er_hist->Fill(this_g1overg2er);
  g1overg2er_hist->GetXaxis()->SetTitle("#sigma_{G1/G2}");
  g1overg2er_hist->Draw();
  sprintf(savename,"plots/M%s/LaserPulsing/Scan%s/g1overg2er_M%s.png", module_str.c_str(), scan_id.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

}