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

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 pede file" << endl;
    cout << "arg 4: specify data file" << endl;
    cout << " " << endl;
    exit(1);
  }

  string module_str = argv[1];
  string data_loc = argv[2];
  string pede_file = argv[3];
  string data_file = argv[4];

  jungfrauFile *thisfile = new jungfrauFile();

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

  jungfrauPedestal *pedestalObject = new jungfrauPedestal();
  pedestalObject->pedestalSetNFrames(1000); // using 1000 frames, rolling window

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

  TH2F* no_g0g1_calib_map = new TH2F("no_g0g1_calib_map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  char savename[128];

  int createHistoFile = 0;
  int createPerStepHistos = 0;

  // open pede file
  sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str());
  thisfile->open((char*)savename, 0);
  
  // calculate pixel mask
  pixelMaskObject->initialisePixelMask(pixel_mask);

  int framecounter = 0;
  while (thisfile->readNextFrame()) {
    framecounter++;
  }
  thisfile->rewind();
  cout << "read " << framecounter << " frames" << endl;

  for (int i = 0; i < 1000; i++) {
    thisfile->readNextFrame();
    pixelMaskObject->maskIfGainNot(0, thisfile->getFrameDataHandle(), pixel_mask);
  }
  for (int i = 0; i < 1000; i++) {
    thisfile->readNextFrame();
    pixelMaskObject->maskIfGainNot(1, thisfile->getFrameDataHandle(), pixel_mask);
  }
  for (int i = 0; i < 999; i++) {
    thisfile->readNextFrame();
    pixelMaskObject->maskIfGainNot(3, thisfile->getFrameDataHandle(), pixel_mask);
  }
  thisfile->rewind();

  sprintf(savename,"plots/M%s/DirectBeam/pixelmask.png", module_str.c_str());
  pixelMaskObject->plotPixelMask(pixel_mask,savename);
  cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;

  // caluclate pedestals
  for (int i = 0; i < 1000; i++) {
    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 i = 0; i < 1000; i++) {
    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 i = 0; i < 999; i++) {
    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();

  thisfile->close();
  
  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);

  TCanvas* c1 = new TCanvas("c1","");

  mapcanvas->cd();

  pedestalsG0->GetXaxis()->SetTitle("Column");
  pedestalsG0->GetYaxis()->SetTitle("Row");
  pedestalsG0->GetYaxis()->SetTitleOffset(0.7);
  pedestalsG0->Draw("colz");
  sprintf(savename,"plots/M%s/DirectBeam/pedeG0_M%s.png", module_str.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/DirectBeam/pedeG1_M%s.png", module_str.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/DirectBeam/pedeG2_M%s.png", module_str.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/DirectBeam/pedeRMSG0_M%s.png", module_str.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/DirectBeam/pedeRMSG1_M%s.png", module_str.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/DirectBeam/pedeRMSG2_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  jungfrauPedestal *pedestalObjectG0 = new jungfrauPedestal();
  pedestalObjectG0->pedestalSetNFrames(1000); // using 1000 frames, rolling window

  jungfrauPedestal *pedestalObjectG1 = new jungfrauPedestal();
  pedestalObjectG1->pedestalSetNFrames(1000); // using 1000 frames, rolling window

  int nsteps = 0;
  if (module_str == "094" || module_str == "066") {
    nsteps = 79;
  } else {
    nsteps = 40;
  }

  // 2D hist for every point, to save 
  // mean, error and number of fills
  // can be loaded from a file
  TH2F* avg_adc_g0_map [nsteps];
  TH2F* avg_adc_g1_map [nsteps];
  TH2F* avg_adcer_g0_map [nsteps];
  TH2F* avg_adcer_g1_map [nsteps];

  // creating the histograms
  if (createHistoFile) {

    for (int j = 0; j < nsteps; j++) {
      sprintf(savename,"avg_adc_g0_map_%d", j);
      avg_adc_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
      sprintf(savename,"avg_adc_g1_map_%d", j);
      avg_adc_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
      
      sprintf(savename,"avg_adcer_g0_map_%d", j);
      avg_adcer_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
      sprintf(savename,"avg_adcer_g1_map_%d", j);
      avg_adcer_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
    }

    int nfiles = ceil(nsteps/10.);

    // total frame number check
    int totalFrames = 0;

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

      sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str());
      thisfile->open((char*)savename, filei);

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

      thisfile->close();
    }

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

    int framecounter = 0;
    for (int filei = 0; filei < nfiles; 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()) {

	if (module_str == "028") {
	  if (framecounter%1000 < 100) { // shutter skip
	    framecounter++;
	    continue;
	  }
	}

	pedestalObjectG0->addG0FrameToPedestalCalculation(thisfile->getFrameDataHandle());
	pedestalObjectG1->addG1FrameToPedestalCalculation(thisfile->getFrameDataHandle());

	framecounter++;

	if (framecounter%1000 == 0) {
	  
	  for (int i = 0; i < NCH; i++) {
	    
	    avg_adc_g0_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG0->pedestalOfChannel(i));
	    avg_adcer_g0_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG0->semOfChannel(i));
	    avg_adc_g1_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG1->pedestalOfChannel(i));
	    avg_adcer_g1_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG1->semOfChannel(i));
	    
	  }
	  
	  pedestalObjectG0->pedestalClear();
	  pedestalObjectG1->pedestalClear();

	  if (createPerStepHistos == 1) {
	    avg_adc_g0_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
	    avg_adc_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
	    avg_adc_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
	    avg_adc_g0_map[framecounter/1000 -1]->Draw("colz");
	    avg_adc_g0_map[framecounter/1000 -1]->SetMinimum(0);
	    sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adc_g0_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
	    mapcanvas->SaveAs((const char *)(savename));
    
	    avg_adc_g1_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
	    avg_adc_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
	    avg_adc_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
	    avg_adc_g1_map[framecounter/1000 -1]->Draw("colz");
	    avg_adc_g1_map[framecounter/1000 -1]->SetMinimum(0);
	    sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adc_g1_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
	    mapcanvas->SaveAs((const char *)(savename));

	    avg_adcer_g0_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
	    avg_adcer_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
	    avg_adcer_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
	    avg_adcer_g0_map[framecounter/1000 -1]->Draw("colz");
	    avg_adcer_g0_map[framecounter/1000 -1]->SetMinimum(0);
	    sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adcer_g0_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
	    mapcanvas->SaveAs((const char *)(savename));
	  
	    avg_adcer_g1_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
	    avg_adcer_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
	    avg_adcer_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
	    avg_adcer_g1_map[framecounter/1000 -1]->Draw("colz");
	    avg_adcer_g1_map[framecounter/1000 -1]->SetMinimum(0);
	    sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adcer_g1_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
	    mapcanvas->SaveAs((const char *)(savename));
	  }

	}

      }

      thisfile->close();
    } // end of files


    // save histograms
    sprintf(savename,"/mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/M%s_CalibAna/DB_histos_M%s.root", module_str.c_str(), module_str.c_str());
    TFile* saved_file = new TFile((const char *)(savename),"RECREATE");

    for (int i = 0; i < nsteps; i++) {
      avg_adc_g0_map[i]->Write();
      avg_adc_g1_map[i]->Write();
      avg_adcer_g0_map[i]->Write();
      avg_adcer_g1_map[i]->Write();
    }
    saved_file->Close();
  } else {
    // load histos
    cout << "LOADING HISTOS" << endl;
    sprintf(savename,"/mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/M%s_CalibAna/DB_histos_M%s.root", module_str.c_str(), module_str.c_str());
    TFile* saved_file = new TFile((const char *)(savename),"READ");

    for (int i = 0; i < nsteps; i++) {
      sprintf(savename,"avg_adc_g0_map_%d", i);
      avg_adc_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
      sprintf(savename,"avg_adc_g1_map_%d", i);
      avg_adc_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename));

      sprintf(savename,"avg_adcer_g0_map_%d", i);
      avg_adcer_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
      sprintf(savename,"avg_adcer_g1_map_%d", i);
      avg_adcer_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
    }
  }

  TH1F* g0overg1hist = new TH1F("g0overg1hist","",100,-34,-24);
  TH2F* g0map = new TH2F("g0map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* g1map = new TH2F("g1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* g0overg1map = new TH2F("g0overg1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* g0ermap = new TH2F("g0ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* g1ermap = new TH2F("g1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH1F* g0overg1erhist = new TH1F("g0overg1erhist","",100,0,2);
  TH2F* g0overg1ermap = new TH2F("g0overg1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  TH1F* g0overg1hist_isEdge = new TH1F("g0overg1hist_isEdge","",100,-34,-24);
  TH1F* g0overg1hist_isInnerEdge = new TH1F("g0overg1hist_isInnerEdge","",100,-34,-24);
  TH1F* g0overg1hist_isDouble = new TH1F("g0overg1hist_isDouble","",100,-34,-24);
  TH1F* g0overg1hist_isNextToDouble = new TH1F("g0overg1hist_isNextToDouble","",100,-34,-24);
  TH1F* g0overg1hist_isQuad = new TH1F("g0overg1hist_isQuad","",100,-34,-24);
  TH1F* g0overg1hist_isBulk = new TH1F("g0overg1hist_isBulk","",100,-34,-24);

  TH1F* range0maxhist = new TH1F("range0maxhist","",100,-10.5,89.5);
  TH1F* range1minhist = new TH1F("range1minhist","",100,-10.5,89.5);
  TH2F* range0maxmap = new TH2F("range0maxmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* range1minmap = new TH2F("range1minmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

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

  for (int i = 0; i < NCH; i++) {

    vector<double> r0_adc;
    vector<double> r0_adcer;
    vector<double> r0_filter;
    vector<double> r0_ferr;

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

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

      int current = 0;
      if (nsteps == 40) {
	current = (j+1)*2;
      } else if (nsteps == 79) {
	current = j+2;
      }

      double this_g0avg = avg_adc_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
      double this_g1avg = avg_adc_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
      double this_g0avger = avg_adcer_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
      double this_g1avger = avg_adcer_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);

      if (this_g0avg != 0) {
	r0_filter.push_back(current);
	r0_ferr.push_back(0);
	r0_adc.push_back(this_g0avg);
	r0_adcer.push_back(this_g0avger);
      }
      if (this_g1avg != 0) {
	r1_filter.push_back(current);
	r1_ferr.push_back(0);
	r1_adc.push_back(this_g1avg);
	r1_adcer.push_back(this_g1avger);
      }
    }

    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;

    // define graphs
    if (r0_adc.size() > 1) {
      grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_ferr[0]),&(r0_adcer[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_adcer[0]));
      grap_g1->SetMarkerStyle(20);
      grap_g1->SetMarkerColor(kGreen+2);
      grap_g1->SetLineColor(kGreen+2);
    }

    // plot the datapoints
    if (r1_adc.size() > 1) {

      if ((i >= 58000 && i < 58000+10) || // bulk
	  (i >= 10 && i < 10+10) || // edge
	  (i >= 1024+10 && i < 1024+10+10) || // inner edge
	  (i >= (256*1024)+10 && i < (256*1024)+10+10) || // double
	  (i >= (257*1024)+10 && i < (257*1024)+10+10) || // next to double
	  (i == (255*1024)+255) // quad
	  ) {

	string pixel_type = "x";
	if (i >= 58000 && i < 58000+10) {
	  pixel_type = "b";
	} else if (i >= 10 && i < 10+10) {
	  pixel_type = "e";
	} else if (i >= 1024+10 && i < 1024+10+10) {
	  pixel_type = "ie";
	} else if (i >= (256*1024)+10 && i < (256*1024)+10+10) {
	  pixel_type = "d";
	} else if (i >= (257*1024)+10 && i < (257*1024)+10+10) {
	  pixel_type = "ntd";
	} else if (i == (255*1024)+255) {
	  pixel_type = "q";
	}

	grap_g1->SetMinimum(0);
	grap_g1->SetMaximum(17000);
	grap_g1->GetXaxis()->SetTitle("Tube current [mA]");
	grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
	grap_g1->GetYaxis()->SetTitleOffset(0.9);
	grap_g1->GetXaxis()->SetLimits(0,85);
	grap_g1->Draw("AP");
	if (r0_adc.size() > 1) {
	  grap_g0->Draw("P");
	}
	sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_M%s.png", module_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) {

      rangemin_g0 = *min_element(r0_filter.begin(),r0_filter.end());
      rangemax_g0 = highestPointBeforeSwitching(r0_filter,r1_filter);
      range0maxhist->Fill(rangemax_g0);
      range0maxmap->Fill(i%NC,i/NC,rangemax_g0);

      if (rangemax_g0 > rangemin_g0) {

	fit_g0 = new TF1("fit_g0","[0]+[1]*x",rangemin_g0,rangemax_g0);
	fit_g0->SetParameter(0, 100.);
	fit_g0->SetParameter(1, 1);
	fit_g0->SetLineColor(kBlue);
	fit_g0->SetParName(0,"G0 const");
	fit_g0->SetParName(1,"G0 grad");
	grap_g0->Fit(fit_g0,"QR","");

	g0map->Fill(i%NC,i/NC,fit_g0->GetParameter(1));
	g0ermap->Fill(i%NC,i/NC,fit_g0->GetParError(1));

	if ((i >= 58000 && i < 58000+10) || // bulk
	    (i >= 10 && i < 10+10) || // edge
	    (i >= 1024+10 && i < 1024+10+10) || // inner edge
	    (i >= (256*1024)+10 && i < (256*1024)+10+10) || // double
	    (i >= (257*1024)+10 && i < (257*1024)+10+10) || // next to double
	    (i == (255*1024)+255) // quad
	    ) {

	  string pixel_type = "x";
	  if (i >= 58000 && i < 58000+10) {
	    pixel_type = "b";
	  } else if (i >= 10 && i < 10+10) {
	    pixel_type = "e";
	  } else if (i >= 1024+10 && i < 1024+10+10) {
	    pixel_type = "ie";
	  } else if (i >= (256*1024)+10 && i < (256*1024)+10+10) {
	    pixel_type = "d";
	  } else if (i >= (257*1024)+10 && i < (257*1024)+10+10) {
	    pixel_type = "ntd";
	  } else if (i == (255*1024)+255) {
	    pixel_type = "q";
	  }

	  grap_g0->SetMinimum(0);
	  grap_g0->SetMaximum(17000);
	  grap_g0->GetXaxis()->SetTitle("Tube current [mA]");
	  grap_g0->GetYaxis()->SetTitle("ADC [ADU]");
	  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");
	  fit_g0->Draw("same");
	  mapcanvas->Update();
	  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);
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));

	  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_adcer[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_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_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_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);

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

	  norm_g0->GetXaxis()->SetTitle("Tube current [mA]");
	  norm_g0->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	  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");
	  lin_g0_p02pc->Draw("same");
	  lin_g0_p05pc->Draw("same");
	  lin_g0_m02pc->Draw("same");
	  lin_g0_m05pc->Draw("same");
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));
	  norm_g0->SetMinimum(-100);
	  norm_g0->SetMaximum(100);
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_zoom_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));

	  delete norm_g0;
	  delete flat_g0;
	}
      }
    }

    if (r1_adc.size() > 1) {

      rangemin_g1 = lowestPointAfterSwitching(r1_filter,r0_filter);
      range1minhist->Fill(rangemin_g1);
      range1minmap->Fill(i%NC,i/NC,rangemin_g1);
      rangemax_g1 = *max_element(r1_filter.begin(),r1_filter.end());

      if (rangemax_g1 > rangemin_g1) {

	fit_g1 = new TF1("fit_g1","[0]+[1]*x",rangemin_g1,rangemax_g1);
	fit_g1->SetParameter(0, 1000.);
	fit_g1->SetParameter(1, -0.1);
	fit_g1->SetLineColor(kGreen+2);
	fit_g1->SetParName(0,"G1 const");
	fit_g1->SetParName(1,"G1 grad");
	grap_g1->Fit(fit_g1,"QR","");

	g1map->Fill(i%NC,i/NC,fit_g1->GetParameter(1));
	g1ermap->Fill(i%NC,i/NC,fit_g1->GetParError(1));

	if ((i >= 58000 && i < 58000+10) || // bulk
	    (i >= 10 && i < 10+10) || // edge
	    (i >= 1024+10 && i < 1024+10+10) || // inner edge
	    (i >= (256*1024)+10 && i < (256*1024)+10+10) || // double
	    (i >= (257*1024)+10 && i < (257*1024)+10+10) || // next to double
	    (i == (255*1024)+255) // quad
	    ) {
	    
	  string pixel_type = "x";
	  if (i >= 58000 && i < 58000+10) {
	    pixel_type = "b";
	  } else if (i >= 10 && i < 10+10) {
	    pixel_type = "e";
	  } else if (i >= 1024+10 && i < 1024+10+10) {
	    pixel_type = "ie";
	  } else if (i >= (256*1024)+10 && i < (256*1024)+10+10) {
	    pixel_type = "d";
	  } else if (i >= (257*1024)+10 && i < (257*1024)+10+10) {
	    pixel_type = "ntd";
	  } else if (i == (255*1024)+255) {
	    pixel_type = "q";
	  }

	  grap_g1->SetMinimum(0);
	  grap_g1->SetMaximum(17000);
	  grap_g1->GetXaxis()->SetTitle("Tube current [mA]");
	  grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
	  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");
	  fit_g1->Draw("same");
	  mapcanvas->Update();
	  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);
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1_M%s.png", module_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]));
	  }

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

	  TF1* flat_g1 = new TF1("flat_g1","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_p05pc = new TF1("lin_g1_p05pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_p05pc->SetParameter(0,fit_g1->GetParameter(0)/200.);
	  lin_g1_p05pc->SetParameter(1,fit_g1->GetParameter(1)/200.);
	  lin_g1_p05pc->SetLineColor(kOrange+1);

	  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_m05pc = new TF1("lin_g1_m05pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
	  lin_g1_m05pc->SetParameter(0,fit_g1->GetParameter(0)/-200.);
	  lin_g1_m05pc->SetParameter(1,fit_g1->GetParameter(1)/-200.);
	  lin_g1_m05pc->SetLineColor(kOrange+1);

	  norm_g1->GetXaxis()->SetTitle("Tube current [mA]");
	  norm_g1->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	  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");
	  lin_g1_p02pc->Draw("same");
	  lin_g1_p05pc->Draw("same");
	  lin_g1_m02pc->Draw("same");
	  lin_g1_m05pc->Draw("same");
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1norm_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));
	  norm_g1->SetMinimum(-50);
	  norm_g1->SetMaximum(50);
	  sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1norm_zoom_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));

	  delete norm_g1;
	  delete flat_g1; 
	}
      }
    }

    // get ratio measurements
    if (rangemax_g0 > rangemin_g0 && rangemax_g1 > rangemin_g1) {

      double this_g0overg1 = fit_g0->GetParameter(1)/fit_g1->GetParameter(1);
      double this_g0overg1er = abs(fit_g0->GetParameter(1)/fit_g1->GetParameter(1))*sqrt(pow((fit_g0->GetParError(1)/fit_g0->GetParameter(1)),2) + pow((fit_g1->GetParError(1)/fit_g1->GetParameter(1)),2));

      g0overg1map->Fill(i%NC,i/NC,this_g0overg1);
      g0overg1hist->Fill(this_g0overg1);
      g0overg1ermap->Fill(i%NC,i/NC,this_g0overg1er);
      g0overg1erhist->Fill(this_g0overg1er);

      if (isEdge(i)) {
	g0overg1hist_isEdge->Fill(this_g0overg1);
      }
      if (isInnerEdge(i)) {
	g0overg1hist_isInnerEdge->Fill(this_g0overg1);
      }
      if (isDouble(i)) {
	g0overg1hist_isDouble->Fill(this_g0overg1);
      }
      if (isNextToDouble(i)) {
	g0overg1hist_isNextToDouble->Fill(this_g0overg1);
      }
      if (isQuad(i)) {
	g0overg1hist_isQuad->Fill(this_g0overg1);
      }
      if (isBulk(i)) {
	g0overg1hist_isBulk->Fill(this_g0overg1);
      }
    } else {
      no_g0g1_calib_map->Fill(i%NC,i/NC,1);
    }

    delete fit_g0;
    delete fit_g1;
    delete grap_g0;
    delete grap_g1;
  }

  mapcanvas->SetLeftMargin(0.1);
  mapcanvas->SetRightMargin(0.13);

  range1minhist->GetXaxis()->SetTitle("Start of gain 1 fit [ADU]");
  range1minhist->Draw();
  sprintf(savename,"plots/M%s/DirectBeam/range1minhist_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  range0maxhist->GetXaxis()->SetTitle("End of gain 0 fit [ADU]");
  range0maxhist->Draw();
  sprintf(savename,"plots/M%s/DirectBeam/range0maxhist_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

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

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

  g1map->GetXaxis()->SetTitle("Column");
  g1map->GetYaxis()->SetTitle("Row");
  g1map->GetYaxis()->SetTitleOffset(0.7);
  g1map->Draw("colz");
  g1map->GetZaxis()->SetRangeUser(-40,0);
  sprintf(savename,"plots/M%s/DirectBeam/g1map_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g0map->GetXaxis()->SetTitle("Column");
  g0map->GetYaxis()->SetTitle("Row");
  g0map->GetYaxis()->SetTitleOffset(0.7);
  g0map->Draw("colz");
  g0map->GetZaxis()->SetRangeUser(200,800);
  sprintf(savename,"plots/M%s/DirectBeam/g0map_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g1ermap->GetXaxis()->SetTitle("Column");
  g1ermap->GetYaxis()->SetTitle("Row");
  g1ermap->GetYaxis()->SetTitleOffset(0.7);
  g1ermap->Draw("colz");
  g1ermap->GetZaxis()->SetRangeUser(0,0.7);
  sprintf(savename,"plots/M%s/DirectBeam/g1ermap_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g0ermap->GetXaxis()->SetTitle("Column");
  g0ermap->GetYaxis()->SetTitle("Row");
  g0ermap->GetYaxis()->SetTitleOffset(0.7);
  g0ermap->Draw("colz");
  g0ermap->GetZaxis()->SetRangeUser(2,8);
  sprintf(savename,"plots/M%s/DirectBeam/g0ermap_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g0overg1map->GetXaxis()->SetTitle("Column");
  g0overg1map->GetYaxis()->SetTitle("Row");
  g0overg1map->GetYaxis()->SetTitleOffset(0.7);
  g0overg1map->Draw("colz");
  g0overg1map->GetZaxis()->SetRangeUser(-35,-20);
  sprintf(savename,"plots/M%s/DirectBeam/g0overg1map_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g0overg1ermap->GetXaxis()->SetTitle("Column");
  g0overg1ermap->GetYaxis()->SetTitle("Row");
  g0overg1ermap->GetYaxis()->SetTitleOffset(0.7);
  g0overg1ermap->Draw("colz");
  g0overg1ermap->GetZaxis()->SetRangeUser(0,1);
  sprintf(savename,"plots/M%s/DirectBeam/g0overg1ermap_M%s.png", module_str.c_str(), module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

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

  c1->cd();

  g0overg1hist->GetXaxis()->SetTitle("G0 / G1");
  g0overg1hist->Draw();
  sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist_M%s.png", module_str.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g0overg1erhist->GetXaxis()->SetTitle("#sigma_{G0 / G1}");
  g0overg1erhist->Draw();
  sprintf(savename,"plots/M%s/DirectBeam/g0overg1erhist_M%s.png", module_str.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g0overg1hist_isEdge->SetLineColor(kBlue);
  g0overg1hist_isInnerEdge->SetLineColor(kCyan);
  g0overg1hist_isDouble->SetLineColor(kGreen+2);
  g0overg1hist_isNextToDouble->SetLineColor(kRed);
  g0overg1hist_isQuad->SetLineColor(kOrange);

  g0overg1hist_isEdge->Scale(1./g0overg1hist_isEdge->GetEntries());
  g0overg1hist_isInnerEdge->Scale(1./g0overg1hist_isInnerEdge->GetEntries());
  g0overg1hist_isDouble->Scale(1./g0overg1hist_isDouble->GetEntries());
  g0overg1hist_isNextToDouble->Scale(1./g0overg1hist_isNextToDouble->GetEntries());
  g0overg1hist_isQuad->Scale(1./g0overg1hist_isQuad->GetEntries());
  g0overg1hist_isBulk->Scale(1./g0overg1hist_isBulk->GetEntries());

  TLegend *leg = new TLegend(0.62,0.6,0.93,0.93);
  leg->AddEntry(g0overg1hist_isBulk, "Normal", "l");
  leg->AddEntry(g0overg1hist_isDouble, "Double", "l");
  leg->AddEntry(g0overg1hist_isNextToDouble, "Next to D", "l");
  leg->AddEntry(g0overg1hist_isEdge, "Edge", "l");
  leg->AddEntry(g0overg1hist_isInnerEdge, "Inner E", "l");

  g0overg1hist_isDouble->GetXaxis()->SetTitle("G0 / G1");
  g0overg1hist_isDouble->GetYaxis()->SetTitle("Normalised");
  g0overg1hist_isDouble->GetYaxis()->SetTitleOffset(1.3);
  g0overg1hist_isDouble->SetMinimum(0.0);
  g0overg1hist_isDouble->SetMaximum(0.1);
  g0overg1hist_isDouble->Draw();
  g0overg1hist_isEdge->Draw("same");
  g0overg1hist_isInnerEdge->Draw("same");
  g0overg1hist_isNextToDouble->Draw("same");
  g0overg1hist_isBulk->Draw("same");
  leg->Draw("same");
  sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist_perType_M%s.png", module_str.c_str(), module_str.c_str());
  c1->SaveAs((const char *)(savename));

  sprintf(savename,"data/M%s/DB_ratio_M%s.root", module_str.c_str(), module_str.c_str());
  TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE");
  g0overg1map->Write();
  g0overg1ermap->Write();
  saved_file2->Close();

}