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

double highestPointBeforeSwitching(const vector<double> &lower_filter, const vector<double> &higher_filter) {

  // find the highest value in lower_filter that is lower than all entries in higher_filter
  double highest_point_before_switching = *min_element(lower_filter.begin(),lower_filter.end())-1.;

  double lowest_entry_in_higher_filter = *min_element(higher_filter.begin(),higher_filter.end());
  for(vector<double>::const_iterator it = lower_filter.begin(); it != lower_filter.end(); ++it) {
    if (*it < lowest_entry_in_higher_filter) {
      if (*it > highest_point_before_switching) {
	highest_point_before_switching = *it;
      }
    }
  }

  return highest_point_before_switching;

}

double lowestPointAfterSwitching(const vector<double> &higher_filter, const vector<double> &lower_filter) {

  // find the lowest value in higher_filter that is higher than all entries in lower_filter
  double lowest_point_after_switching = *max_element(higher_filter.begin(),higher_filter.end())+1;

  double highest_entry_in_lower_filter = *max_element(lower_filter.begin(),lower_filter.end());
  for(vector<double>::const_iterator it = higher_filter.begin(); it != higher_filter.end(); ++it) {
    if (*it > highest_entry_in_lower_filter) {
      if (*it < lowest_point_after_switching) {
	lowest_point_after_switching = *it;
      }
    }
  }

  return lowest_point_after_switching;

}

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

  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);
  if (module_str == "028") {
    pixelMaskObject->maskChip(1, pixel_mask);
    pixelMaskObject->maskChip(2, pixel_mask);
    pixelMaskObject->maskChip(5, pixel_mask);
    pixelMaskObject->maskChip(6, pixel_mask);
  } else if (module_str == "006") {
    pixelMaskObject->maskChip(6, pixel_mask);
    pixelMaskObject->maskSupercolumn(4,1, pixel_mask);
    pixelMaskObject->maskSupercolumn(4,2, pixel_mask);
  } else if (module_str == "094") {
    pixelMaskObject->maskChip(6,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);

  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 = 80;
  } 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 = nsteps/10;
    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++) {
	    if (pixel_mask[i] == true) {
	      
	      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);

  for (int i = 0; i < NCH; i++) {
    if (pixel_mask[i] == true) {

      vector<double> r0_adc;
      vector<double> r0_adcer;
      vector<double> r0_filter;
      vector<double> r0_filterer;
      vector<double> r1_adc;
      vector<double> r1_adcer;
      vector<double> r1_filter;
      vector<double> r1_filterer;

      for (int j = 0; j < nsteps; j++) {
	int current = 0;
	if (nsteps == 40) {
	  current = (j+1)*2;
	} else if (nsteps == 80) {
	  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 && this_g0avger > 0) {
	  r0_filter.push_back(current);
	  r0_filterer.push_back(0);
	  r0_adc.push_back(this_g0avg);
	  r0_adcer.push_back(this_g0avger);
	}
	if (this_g1avg != 0 && this_g1avger > 0) {
	  r1_filter.push_back(current);
	  r1_filterer.push_back(0);
	  r1_adc.push_back(this_g1avg);
	  r1_adcer.push_back(this_g1avger);
	}

      }

      if (r0_adc.size() > 1 && r1_adc.size() > 1) {

	TGraphErrors *grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_filterer[0]),&(r0_adcer[0]));
	TGraphErrors *grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_filterer[0]),&(r1_adcer[0]));

	double rangemin0 = 0;
	if (module_str == "028" || module_str == "021" || module_str == "022" || module_str == "008" || module_str == "044" || module_str == "066" || module_str == "094") {
	  rangemin0 = 2;
	} else if (module_str == "006") {
	  rangemin0 = 4;
	}

	double rangemax0 = highestPointBeforeSwitching(r0_filter,r1_filter);
	range0maxhist->Fill(rangemax0);
	range0maxmap->Fill(i%NC,i/NC,rangemax0);

	TF1 *fit1 = new TF1("fit1","[0]+[1]*x",rangemin0,rangemax0);
	TF1 *fit2 = 0;
	TF1 *fit1_e = 0;
	TF1 *fit2_e = 0;
	
	fit1->SetParameter(0, 100.);
	fit1->SetParameter(1, 1);
	fit1->SetParName(0,"G0 const");
	fit1->SetParName(1,"G0 grad");
	fit1->SetLineColor(kBlue);
	grap_g0->Fit(fit1,"QR","");

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

	fit1_e = new TF1("fit1_e","[0]+[1]*x",0,40);
	fit1_e->SetParameters(fit1->GetParameter(0),fit1->GetParameter(1));
	fit1_e->SetLineColor(kBlue);
	fit1_e->SetLineStyle(2);
      
	double rangemin1 = lowestPointAfterSwitching(r1_filter,r0_filter);
	range1minhist->Fill(rangemin1);
	range1minmap->Fill(i%NC,i/NC,rangemin1);

	double rangemax1 = 80;

	fit2 = new TF1("fit2","[0]+[1]*x",rangemin1,rangemax1);

	fit2->SetParameter(0, 1000.);
	fit2->SetParameter(1, -0.1);
	fit2->SetParName(0,"G1 const");
	fit2->SetParName(1,"G1 grad");
	fit2->SetLineColor(kGreen+2);
	grap_g1->Fit(fit2,"QR","");

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

	fit2_e = new TF1("fit2_e","[0]+[1]*x",0,80);
	fit2_e->SetParameters(fit2->GetParameter(0),fit2->GetParameter(1));
	fit2_e->SetLineColor(kGreen+2);
	fit2_e->SetLineStyle(2);
      
	g0overg1map->Fill(i%NC,i/NC,fit1->GetParameter(1)/fit2->GetParameter(1));
	g0overg1hist->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	g0overg1ermap->Fill(i%NC,i/NC,abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));
	g0overg1erhist->Fill(abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));

	if (isEdge(i)) {
	  g0overg1hist_isEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	}
	if (isInnerEdge(i)) {
	  g0overg1hist_isInnerEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	}
	if (isDouble(i)) {
	  g0overg1hist_isDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	}
	if (isNextToDouble(i)) {
	  g0overg1hist_isNextToDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	}
	if (isQuad(i)) {
	  g0overg1hist_isQuad->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
	}
	if (isBulk(i)) {
	  g0overg1hist_isBulk->Fill(fit1->GetParameter(1)/fit2->GetParameter(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";
	  }

	  mapcanvas->SetLeftMargin(0.13);
	  mapcanvas->SetRightMargin(0.05);
	    
	  grap_g0->SetMarkerStyle(20);
	  grap_g0->SetMarkerColor(kBlue);
	  grap_g0->SetLineColor(kBlue);
	  grap_g1->SetMarkerStyle(20);
	  grap_g1->SetMarkerColor(kGreen+2);
	  grap_g1->SetLineColor(kGreen+2);

	  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->SetMinimum(0);
	  grap_g1->SetMaximum(17000);
	  grap_g1->Draw("AP");
	  grap_g0->Draw("P");
	  fit1->Draw("same");
	  fit1_e->Draw("same");
	  fit2->Draw("same");
	  fit2_e->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);

	  TPaveStats *st1 = (TPaveStats*)grap_g1->FindObject("stats");
	  st1->SetX1NDC(0.6);
	  st1->SetX2NDC(0.94);
	  st1->SetY1NDC(0.18);
	  st1->SetY2NDC(0.37);
	  st1->SetBorderSize(0);
	  st1->SetTextColor(kGreen+2);

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

	  // normalise
	  vector<double> r0_adc_norm;
	  vector<double> r1_adc_norm;
	  for (size_t j = 0; j < r0_adc.size(); j++) {
	    r0_adc_norm.push_back(r0_adc[j] - fit1->Eval(r0_filter[j]));
	  }
	  for (size_t j = 0; j < r1_adc.size(); j++) {
	    r1_adc_norm.push_back(r1_adc[j] - fit2->Eval(r1_filter[j]));
	  }
	  TGraphErrors *grap0_norm = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc_norm[0]),&(r0_filterer[0]),&(r0_adcer[0]));
	  TGraphErrors *grap1_norm = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_filterer[0]),&(r1_adcer[0]));

	  TF1* flat_g0 = new TF1("flat_g0","0",rangemin0, rangemax0);
	  TF1* flat_g1 = new TF1("flat_g1","0",rangemin1, rangemax1);

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

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

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

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

	  TF1* lin_g1_p02pc = new TF1("lin_g1_p02pc","[0]+[1]*x",rangemin1,rangemax1);
	  lin_g1_p02pc->SetParameter(0,fit2->GetParameter(0)/500.);
	  lin_g1_p02pc->SetParameter(1,fit2->GetParameter(1)/500.);
	  lin_g1_p02pc->SetLineColor(kOrange);

	  TF1* lin_g1_p05pc = new TF1("lin_g1_p05pc","[0]+[1]*x",rangemin1,rangemax1);
	  lin_g1_p05pc->SetParameter(0,fit2->GetParameter(0)/200.);
	  lin_g1_p05pc->SetParameter(1,fit2->GetParameter(1)/200.);
	  lin_g1_p05pc->SetLineColor(kOrange+1);

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

	  TF1* lin_g1_m05pc = new TF1("lin_g1_m05pc","[0]+[1]*x",rangemin1,rangemax1);
	  lin_g1_m05pc->SetParameter(0,fit2->GetParameter(0)/-200.);
	  lin_g1_m05pc->SetParameter(1,fit2->GetParameter(1)/-200.);
	  lin_g1_m05pc->SetLineColor(kOrange+1);

	  grap0_norm->SetMarkerColor(kBlue);
	  grap0_norm->SetLineColor(kBlue);
	  flat_g0->SetLineColor(kBlue);
	  grap0_norm->GetXaxis()->SetTitle("Tube current [mA]");
	  grap0_norm->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	  grap0_norm->GetYaxis()->SetTitleOffset(0.9);
	  grap0_norm->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));
	  grap0_norm->SetMinimum(-100);
	  grap0_norm->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));

	  grap1_norm->SetMarkerColor(kGreen+2);
	  grap1_norm->SetLineColor(kGreen+2);
	  flat_g1->SetLineColor(kGreen+2);
	  grap1_norm->GetXaxis()->SetTitle("Tube current [mA]");
	  grap1_norm->GetYaxis()->SetTitle("Normalised ADC [ADU]");
	  grap1_norm->GetYaxis()->SetTitleOffset(0.9);
	  grap1_norm->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));
	  grap1_norm->SetMinimum(-50);
	  grap1_norm->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));

	  mapcanvas->SetLeftMargin(0.1);
	  mapcanvas->SetRightMargin(0.13);
	}
	
	delete fit1;
	delete fit1_e;
	delete fit2;
	delete fit2_e;

      }	
    }
  }

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

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

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

}