// file to take current source scan files, analyse and plot them.

#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonHeader.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonFunctions.h"

#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauFile.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPedestal.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPixelMask.C"

#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"
#include "TLegend.h"

#include <algorithm>
#include <sstream>

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

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

  if (argc != 4) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify data location" << endl;
    cout << "arg 3: specify data file" << endl;
    cout << "eg: ./JFMC_CurrentSourceScan 006 ..." << endl;
    cout << " " << endl;
    exit(1);
  }

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

  jungfrauFile *thisfile = new jungfrauFile();
  char savename[128];
  sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str());
  thisfile->open((char*)savename, 0);

  TH2F *gainmap_all = new TH2F("gainmap_all","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *adcmap_all = new TH2F("adcmap_all","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *gainmap_sel = new TH2F("gainmap_sel","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *adcmap_sel = new TH2F("adcmap_sel","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  TH1F *gain_avg = new TH1F("gain_avg","",100,-1,4);
  TH2F *gainmap_avg = new TH2F("gainmap_avg","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH1F *adc_avg_g1 = new TH1F("adc_avg_g1","",100,0,17000);
  TH2F *adcmap_avg_g1 = new TH2F("adcmap_avg_g1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH1F *adc_avg_g2 = new TH1F("adc_avg_g2","",100,0,17000);
  TH2F *adcmap_avg_g2 = new TH2F("adcmap_avg_g2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

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

  jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
  static int pixel_mask [NCH];
  pixelMaskObject->initialisePixelMask(pixel_mask);
  if (module_str == "032") {
    pixelMaskObject->maskChip(2,pixel_mask);
    pixelMaskObject->maskChip(3,pixel_mask);
    pixelMaskObject->maskChip(4,pixel_mask);
  } else if (module_str == "006") {
    pixelMaskObject->maskChip(6, pixel_mask);
    pixelMaskObject->maskSupercolumn(4,1, pixel_mask);
    pixelMaskObject->maskSupercolumn(4,2, pixel_mask);
  }

  int frame_counter = 0;

  static TH1I *adc_histos_g1 [NCH];
  static TH1I *adc_histos_g2 [NCH];
  static TH1I *gain_histos [NCH];

  for(int i = 0; i < NCH; i++) {
    ostringstream histogramNameStreamg1;
    histogramNameStreamg1 << "adc_histos_g1_" << i;
    adc_histos_g1[i] = new TH1I(histogramNameStreamg1.str().c_str(),"",100,0,17000);
    ostringstream histogramNameStreamg2;
    histogramNameStreamg2 << "adc_histos_g2_" << i;
    adc_histos_g2[i] = new TH1I(histogramNameStreamg2.str().c_str(),"",100,0,17000);
    ostringstream histogramNameStream2;
    histogramNameStream2 << "gain_histos_" << i;
    gain_histos[i] = new TH1I(histogramNameStream2.str().c_str(),"",100,-1,4);
  }

  double filter[35];
  if (module_str == "044") {
    for (int i = 0; i < 5; i++) {
      filter[i] = 0.05+(i*0.05);
    }
    for (int i = 0; i < 9; i++) {
      filter[i+5] = 0.5+(i*0.25);
    }
    for (int i = 0; i < 17; i++) {
      filter[i+5+9] = 5.+(i*2.5);
    }
  } else {
    for (int i = 0; i < 9; i++) {
      filter[i] = 0.05+(i*0.025);
    }
    for (int i = 0; i < 9; i++) {
      filter[i+9] = 0.5+(i*0.25);
    }
    for (int i = 0; i < 17; i++) {
      filter[i+9+9] = 5.+(i*2.5);
    }
  }

  static double adcmeans_g1[NCH][35];
  static double adcmeans_g2[NCH][35];
  static double adcmeanerrs_g1[NCH][35];
  static double adcmeanerrs_g2[NCH][35];

  float pede_ene[1];
  pede_ene[0] = 0.01;

  float pede_val1[NCH];
  float pede_val2[NCH];

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

  for (int j = 0; j < 640; j++) {
    frame_counter++;
    thisfile->readNextFrame();
    pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
  }
  cout << pedestalObject->pedestalOfChannel(0) << endl;
  cout << pedestalObject->rmsOfChannel(0) << endl;
  pedestalObject->pedestalClear();

  for (int j = 0; j < 640; j++) {
    frame_counter++;
    thisfile->readNextFrame();
    pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
  }
  for (int i = 0; i < NCH; i++) {
    pede_val1[i] = pedestalObject->pedestalOfChannel(i);
  }
  cout << pedestalObject->pedestalOfChannel(0) << endl;
  cout << pedestalObject->rmsOfChannel(0) << endl;
  pedestalObject->pedestalClear();

  for (int j = 0; j < 640; j++) {
    frame_counter++;
    thisfile->readNextFrame();
    pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
  }
  for (int i = 0; i < NCH; i++) {
    pede_val2[i] = pedestalObject->pedestalOfChannel(i);
  }
  cout << pedestalObject->pedestalOfChannel(0) << endl;
  cout << pedestalObject->rmsOfChannel(0) << endl;
  pedestalObject->pedestalClear();

  int j = 0;
  int f = 1;
  int f_max = 0;
  if (module_str == "032" || module_str == "021" || module_str == "022") {
    f_max = 3;
  } else if (module_str == "006" || module_str == "008" || module_str == "044") {
    f_max = 4;
  }
  
  while (f < f_max) {
    while (thisfile->readNextFrame()) {
      frame_counter++;
      j++;

      uint16_t* imagedptr = thisfile->getFrameDataHandle();
      for (int i = 0; i < NCH; i++) {

	if ((j-1) < 10 || (j > 8075 && j < 8085)) {
	  if (pixel_mask[i] == 1) {
	    uint16_t adc = imagedptr[i]&0x3fff;
	    uint16_t gain = (imagedptr[i]&0xc000) >> 14;
	    adcmap_all->Fill(i%NC,i/NC,adc);
	    gainmap_all->Fill(i%NC,i/NC,gain);
	  }
	}

	if (((i/NC <= 255) && (i%64 == (j-1)%64)) ||
	    ((i/NC >= 256) &&((-1*(i-524287))%64 == (j-1)%64))) {

	  if (pixel_mask[i] == 1) {
	    uint16_t adc = imagedptr[i]&0x3fff;
	    uint16_t gain = (imagedptr[i]&0xc000) >> 14;

	    gain_histos[i]->Fill(gain);
	    if (gain == 1) {
	      adc_histos_g1[i]->Fill(adc);
	    } else if (gain == 3) {
	      adc_histos_g2[i]->Fill(adc);
	    }

	    if ((j-1) < 10 || (j > 8075 && j < 8085)) {
	      adcmap_sel->Fill(i%NC,i/NC,adc);
	      gainmap_sel->Fill(i%NC,i/NC,gain);
	    }

	  }
	}
      }
      if ((j-1) < 10 || (j > 8075 && j < 8085)) {
	mapcanvas->cd();
	adcmap_sel->GetXaxis()->SetTitle("Column");
	adcmap_sel->GetYaxis()->SetTitle("Row");
	adcmap_sel->GetYaxis()->SetTitleOffset(0.7);
	adcmap_sel->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/adcmap_%i_sel.png", module_str.c_str(), j-1);
	mapcanvas->SaveAs((const char *)(savename));
	adcmap_all->GetXaxis()->SetTitle("Column");
	adcmap_all->GetYaxis()->SetTitle("Row");
	adcmap_all->GetYaxis()->SetTitleOffset(0.7);
	adcmap_all->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/adcmap_%i_all.png", module_str.c_str(), j-1);
	mapcanvas->SaveAs((const char *)(savename));
      
	gainmap_sel->GetXaxis()->SetTitle("Column");
	gainmap_sel->GetYaxis()->SetTitle("Row");
	gainmap_sel->GetYaxis()->SetTitleOffset(0.7);
	gainmap_sel->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/gainmap_%i_sel.png", module_str.c_str(), j-1);
	mapcanvas->SaveAs((const char *)(savename));
	gainmap_all->GetXaxis()->SetTitle("Column");
	gainmap_all->GetYaxis()->SetTitle("Row");
	gainmap_all->GetYaxis()->SetTitleOffset(0.7);
	gainmap_all->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/gainmap_%i_all.png", module_str.c_str(), j-1);
	mapcanvas->SaveAs((const char *)(savename));
      
	adcmap_sel->Reset();
	adcmap_all->Reset();
	gainmap_sel->Reset();
	gainmap_all->Reset();
      }

      // catch the last frame (bad ctrl-c exit)
      int last_frame = 0;
      if (module_str == "032" || module_str == "021" || module_str == "022") {
	last_frame = 19199;
      } else if (module_str == "006" || module_str == "008") {
	last_frame = 24319;
      } else if (module_str == "044") {
	last_frame = 24319;
      }

      if (j%640 == 0 || j == last_frame) {
	if (j == last_frame) {
	  j = last_frame+1;
	}
	
	for (int i = 0; i < NCH; i++) {
	  if (pixel_mask[i] == 1) {
	    adcmeans_g1[i][j/640-1] = adc_histos_g1[i]->GetMean();
	    adcmeans_g2[i][j/640-1] = adc_histos_g2[i]->GetMean();
	    adcmeanerrs_g1[i][j/640-1] = adc_histos_g1[i]->GetMeanError();
	    adcmeanerrs_g2[i][j/640-1] = adc_histos_g2[i]->GetMeanError();
	    
	    if (gain_histos[i]->GetEntries() > 0) {
	      gain_avg->Fill(gain_histos[i]->GetMean());
	      gainmap_avg->Fill(i%NC,i/NC,gain_histos[i]->GetMean());
	    }
	    gain_histos[i]->Reset();
	    if (adc_histos_g1[i]->GetEntries() > 0) {
	      adc_avg_g1->Fill(adc_histos_g1[i]->GetMean());
	      adcmap_avg_g1->Fill(i%NC,i/NC,adc_histos_g1[i]->GetMean());
	    }
	    adc_histos_g1[i]->Reset();
	    if (adc_histos_g2[i]->GetEntries() > 0) {
	      adc_avg_g2->Fill(adc_histos_g2[i]->GetMean());
	      adcmap_avg_g2->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMean());
	    }
	    adc_histos_g2[i]->Reset();
	  }
	}
	c1->cd();
	gain_avg->GetXaxis()->SetTitle("Average gain");
	gain_avg->GetXaxis()->SetTitleOffset(1.1);
	gain_avg->GetYaxis()->SetTitleOffset(1.5);
	gain_avg->Draw();
	sprintf(savename,"plots/M%s/CurrentSource/perStep/gain_avg_%i.png", module_str.c_str(), j/640);
	c1->SaveAs((const char *)(savename));
	mapcanvas->cd();
	gainmap_avg->GetXaxis()->SetTitle("Column");
	gainmap_avg->GetYaxis()->SetTitle("Row");
	gainmap_avg->GetYaxis()->SetTitleOffset(0.7);
	gainmap_avg->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/perStep/gainmap_avg_%i.png", module_str.c_str(), j/640);
	mapcanvas->SaveAs((const char *)(savename));

	c1->cd();
	adc_avg_g1->GetXaxis()->SetTitle("Average G1 ADC");
	adc_avg_g1->GetXaxis()->SetTitleOffset(1.1);
	adc_avg_g1->GetYaxis()->SetTitleOffset(1.5);
	adc_avg_g1->Draw();
	sprintf(savename,"plots/M%s/CurrentSource/perStep/adc_avg_g1_%i.png", module_str.c_str(), j/640);
	c1->SaveAs((const char *)(savename));
	mapcanvas->cd();
	adcmap_avg_g1->GetXaxis()->SetTitle("Column");
	adcmap_avg_g1->GetYaxis()->SetTitle("Row");
	adcmap_avg_g1->GetYaxis()->SetTitleOffset(0.7);
	adcmap_avg_g1->Draw("colz");
	mapcanvas->Update();
	sprintf(savename,"plots/M%s/CurrentSource/perStep/adcmap_avg_g1_%i.png", module_str.c_str(), j/640);
	mapcanvas->SaveAs((const char *)(savename));

	c1->cd();
	adc_avg_g2->GetXaxis()->SetTitle("Average G2 ADC");
	adc_avg_g2->GetXaxis()->SetTitleOffset(1.1);
	adc_avg_g2->GetYaxis()->SetTitleOffset(1.5);
	adc_avg_g2->Draw();
	sprintf(savename,"plots/M%s/CurrentSource/perStep/adc_avg_g2_%i.png", module_str.c_str(), j/640);
	c1->SaveAs((const char *)(savename));
	mapcanvas->cd();
	adcmap_avg_g2->GetXaxis()->SetTitle("Column");
	adcmap_avg_g2->GetYaxis()->SetTitle("Row");
	adcmap_avg_g2->GetYaxis()->SetTitleOffset(0.7);
	adcmap_avg_g2->Draw("colz");
	sprintf(savename,"plots/M%s/CurrentSource/perStep/adcmap_avg_g2_%i.png", module_str.c_str(), j/640);
	mapcanvas->SaveAs((const char *)(savename));
	gain_avg->Reset();
	gainmap_avg->Reset();
	adc_avg_g1->Reset();
	adcmap_avg_g1->Reset();
	adc_avg_g2->Reset();
	adcmap_avg_g2->Reset();
      }
    
    } // end of while

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

  } // end of file while
  thisfile->close();

  for (int i = 0; i < NCH; i++) {
    delete gain_histos[i];
    delete adc_histos_g1[i];
    delete adc_histos_g2[i];
  }
 
  cout << "frame_counter " << frame_counter << endl;

  TH2F *g1map = new TH2F("g1map","",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);

  TH2F *g2map = new TH2F("g2map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *g2ermap = new TH2F("g2ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  
  TH2F *g1overg2map = new TH2F("g1overg2map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *g1overg2ermap = new TH2F("g1overg2ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  
  TH2F *diffFromLinG1map = new TH2F("diffFromLinG1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F *diffFromLinG2map = new TH2F("diffFromLinG2map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  
  TH1F *g1hist;
  if (module_str == "032") {
    g1hist = new TH1F("g1hist","",100,-12000,-4000);
  } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008" || module_str == "044") {
    g1hist = new TH1F("g1hist","",100,-3000,0);
  }
  TH1F *g1erhist = new TH1F("g1erhist","",100,0,50);

  TH1F *g2hist;
  if (module_str == "032") {
    g2hist = new TH1F("g2hist","",100,-1000,-200);
  } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008" || module_str == "044") {
    g2hist = new TH1F("g2hist","",100,-300,0);
  }
  TH1F *g2erhist = new TH1F("g2erhist","",100,0,10);

  TH1F *g1overg2hist = new TH1F("g1overg2hist","",100,11,15);
  TH1F *g1overg2erhist = new TH1F("g1overg2erhist","",100,0,1);

  TH1F *diffFromLinG1hist = new TH1F("diffFromLinG1hist","",100,-0.05,0.05);
  TH1F *diffFromLinG2hist = new TH1F("diffFromLinG2hist","",100,-0.1,0.1);

  TH1F *g1overg2_isEdge = new TH1F("g1overg2_isEdge","",100,11,15);
  TH1F *g1overg2_isInnerEdge = new TH1F("g1overg2_isInnerEdge","",100,11,15);
  TH1F *g1overg2_isDouble = new TH1F("g1overg2_isDouble","",100,11,15);
  TH1F *g1overg2_isNextToDouble = new TH1F("g1overg2_isNextToDouble","",100,11,15);
  TH1F *g1overg2_isQuad = new TH1F("g1overg2_isQuad","",100,11,15);
  TH1F *g1overg2_isBulk = new TH1F("g1overg2_isBulk","",100,11,15);

  for (int i = 0; i < NCH; i++) {
    if (pixel_mask[i] == 1) {
      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> r3_adc;
      vector<double> r3_filter;
      vector<double> r3_adcerr;
      vector<double> r3_ferr;

      // filter out points at zero
      int m_max = 0;
      if (module_str == "032" || module_str == "021" || module_str == "022" || module_str == "044") {
	m_max = 27;
      } else if (module_str == "006" || module_str == "008") {
	m_max = 35;
      }

      for (int m = 0; m < m_max; m++) {
	if (adcmeans_g1[i][m] != 0) {
	  r1_adc.push_back(adcmeans_g1[i][m]);
	  r1_filter.push_back(filter[m]);
	  r1_adcerr.push_back(adcmeanerrs_g1[i][m]);
	  r1_ferr.push_back(0.);
	}
	if (adcmeans_g2[i][m] != 0) {
	  r3_adc.push_back(adcmeans_g2[i][m]);
	  r3_filter.push_back(filter[m]);
	  r3_adcerr.push_back(adcmeanerrs_g2[i][m]);
	  r3_ferr.push_back(0.);
	}
      }

      float this_pede_val1[1] = {pede_val1[i]};
      float this_pede_val2[1] = {pede_val2[i]};

      TGraph *grap_p1 = new TGraph(1,pede_ene,this_pede_val1);
      TGraph *grap_p2 = new TGraph(1,pede_ene,this_pede_val2);
      
      TGraphErrors *grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_ferr[0]),&(r1_adcerr[0]));
      TGraphErrors *grap_g2 = new TGraphErrors(r3_adc.size(),&(r3_filter[0]),&(r3_adc[0]),&(r3_ferr[0]),&(r3_adcerr[0]));

      TF1 *fit2;
      TF1 *fit3;

      if (r1_adc.size() > 1) {
	double rangemin2 = 0;
	double rangemax2 = 0;
	if (module_str == "032") {
	  rangemin2 = 0.75;
	  rangemax2 = *max_element(r1_filter.begin(),r1_filter.end());
	} else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008" || module_str == "044") {
	  rangemin2 = 0.075;
	  rangemax2 = 2.5;
	}
	fit2 = new TF1("fit2","[0]+[1]*x",rangemin2,rangemax2);
	fit2->SetParameter(0, 10000.);
	fit2->SetParameter(1, -0.1);
	fit2->SetLineColor(kGreen+2);
	fit2->SetParName(0,"G1 const");
	fit2->SetParName(1,"G1 grad");
	grap_g1->Fit(fit2,"QR","");

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

	double g1residsum = 0.;
	double g1residn = 0.;
	for (size_t it = 0; it < r1_adc.size(); it++) {
	  if (r1_filter[it] >= rangemin2) {
	    diffFromLinG1hist->Fill((r1_adc[it] - fit2->Eval(r1_filter[it])) / r1_adc[it]);
	    g1residsum = g1residsum + ((r1_adc[it] - fit2->Eval(r1_filter[it])) / r1_adc[it]);
	    g1residn = g1residn + 1.;
	  }
	}
	diffFromLinG1map->Fill(i%NC,i/NC,g1residsum/g1residn);
      }

      if (r3_adc.size() > 1) {
	double rangemin3 = 0;
	double rangemax3 = 0;
	if (module_str == "032" || module_str == "022" || module_str == "021") {
	  rangemin3 = 5;
	  rangemax3 = *max_element(r3_filter.begin(),r3_filter.end());
	} else if (module_str == "006" || module_str == "008" || module_str == "044") {
	  rangemin3 = 7.5;
	  rangemax3 = *max_element(r3_filter.begin(),r3_filter.end());
	}
	fit3 = new TF1("fit3","[0]+[1]*x",rangemin3, rangemax3);
	fit3->SetParameter(0, 10000.);
	fit3->SetParameter(1, -0.01);
	fit3->SetLineColor(kRed);
	fit3->SetParName(0,"G2 const");
	fit3->SetParName(1,"G2 grad");
	grap_g2->Fit(fit3,"QR+","");

	g2hist->Fill(fit3->GetParameter(1));
	g2map->Fill(i%NC,i/NC,fit3->GetParameter(1));
	g2erhist->Fill(fit3->GetParError(1));
	g2ermap->Fill(i%NC,i/NC,fit3->GetParError(1));

	double g2residsum = 0.;
	double g2residn = 0.;
	for (size_t it = 0; it < r3_adc.size(); it++) {
	  diffFromLinG2hist->Fill((r3_adc[it] - fit3->Eval(r3_filter[it])) / r3_adc[it]);
	  g2residsum = g2residsum + ((r3_adc[it] - fit3->Eval(r3_filter[it])) / r3_adc[it]);
	  g2residn = g2residn + 1.;
	}
	diffFromLinG2map->Fill(i%NC,i/NC,g2residsum/g2residn);
      }

      if (r1_adc.size() > 1 && r3_adc.size() > 1) {
	g1overg2hist->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	g1overg2erhist->Fill(abs(fit2->GetParameter(1)/fit3->GetParameter(1))*sqrt(pow((fit2->GetParError(1)/fit2->GetParameter(1)),2) + pow((fit3->GetParError(1)/fit3->GetParameter(1)),2)));
	g1overg2map->Fill(i%NC,i/NC,fit2->GetParameter(1) / fit3->GetParameter(1));
	g1overg2ermap->Fill(i%NC,i/NC,abs(fit2->GetParameter(1)/fit3->GetParameter(1))*sqrt(pow((fit2->GetParError(1)/fit2->GetParameter(1)),2) + pow((fit3->GetParError(1)/fit3->GetParameter(1)),2)));

	if (isEdge(i)) {
	  g1overg2_isEdge->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	}
	if (isInnerEdge(i)) {
	  g1overg2_isInnerEdge->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	}
	if (isDouble(i)) {
	  g1overg2_isDouble->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	}
	if (isNextToDouble(i)) {
	  g1overg2_isNextToDouble->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	}
	if (isQuad(i)) {
	  g1overg2_isQuad->Fill(fit2->GetParameter(1) / fit3->GetParameter(1));
	}
	if (isBulk(i)) {
	  g1overg2_isBulk->Fill(fit2->GetParameter(1) / fit3->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_p1->SetMarkerColor(kGreen+2);
	grap_p2->SetMarkerColor(kRed);
	grap_p1->SetMarkerStyle(20);
	grap_p2->SetMarkerStyle(20);
	grap_g1->SetMarkerStyle(20);
	grap_g2->SetMarkerStyle(20);
	grap_g1->SetMarkerColor(kGreen+2);
	grap_g2->SetMarkerColor(kRed);
	grap_g1->SetLineColor(kGreen+2);
	grap_g2->SetLineColor(kRed);

	TF1 *fit2_e = new TF1("fit2_e","[0]+[1]*x",0.009,1.1E3);
	TF1 *fit3_e = new TF1("fit3_e","[0]+[1]*x",0.009,1.1E3);

	fit2_e->SetParameters(fit2->GetParameter(0),fit2->GetParameter(1));
	fit2_e->SetLineColor(kGreen+2);
	fit2_e->SetLineStyle(2);
	fit3_e->SetParameters(fit3->GetParameter(0),fit3->GetParameter(1));
	fit3_e->SetLineColor(kRed);
	fit3_e->SetLineStyle(2);

	grap_g1->GetXaxis()->SetLimits(0.009,1.1E2);
	grap_g1->SetMinimum(0);
	grap_g1->SetMaximum(20000);
	grap_g1->GetXaxis()->SetTitle("Integration time [#mus]");
	grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
	grap_g1->GetYaxis()->SetTitleOffset(0.9);
	grap_g1->Draw("AP");
	fit2->Draw("same");
	mapcanvas->Update();
	TPaveStats *st0 = (TPaveStats*)grap_g1->FindObject("stats");
	st0->SetX1NDC(0.2);
	st0->SetX2NDC(0.54);
	st0->SetY1NDC(0.71);
	st0->SetY2NDC(0.90);
	st0->SetBorderSize(0);
	st0->SetTextColor(kGreen+2);
	grap_g2->Draw("P");
	fit3->Draw("same");
	mapcanvas->Update();
	TPaveStats *st1 = (TPaveStats*)grap_g2->FindObject("stats");
	st1->SetX1NDC(0.6);
	st1->SetX2NDC(0.94);
	st1->SetY1NDC(0.71);
	st1->SetY2NDC(0.90);
	st1->SetBorderSize(0);
	st1->SetTextColor(kRed);

	fit2_e->Draw("same");
	fit3_e->Draw("same");
	grap_p1->Draw("P");
	grap_p2->Draw("P");
	mapcanvas->SetLogx();
	sprintf(savename,"plots/M%s/CurrentSource/pixel_%s_%d.png", module_str.c_str(), pixel_type.c_str(), i);
	mapcanvas->SaveAs((const char *)(savename));
	mapcanvas->SetLogx(0);

	mapcanvas->SetLeftMargin(0.1);
	mapcanvas->SetRightMargin(0.13);
	delete fit2_e;
	delete fit3_e;
      }

      delete grap_g1;
      delete grap_g2;
      if (r1_adc.size() > 1) {
	delete fit2;
      }
      if (r3_adc.size() > 1) {
	delete fit3;
      }
    }
  }

  c1->cd();

  g1hist->GetXaxis()->SetTitle("G1 gradient");
  g1hist->GetXaxis()->SetTitleOffset(1.1);
  g1hist->GetYaxis()->SetTitleOffset(1.5);
  g1hist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g1hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g2hist->GetXaxis()->SetTitle("G2 gradient");
  g2hist->GetXaxis()->SetTitleOffset(1.1);
  g2hist->GetYaxis()->SetTitleOffset(1.5);
  g2hist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g2hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g1overg2hist->GetXaxis()->SetTitle("G1 / G2");
  g1overg2hist->GetXaxis()->SetTitleOffset(1.1);
  g1overg2hist->GetYaxis()->SetTitleOffset(1.5);
  g1overg2hist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g1overg2hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g1overg2erhist->GetXaxis()->SetTitle("G1 / G2 uncert");
  g1overg2erhist->GetXaxis()->SetTitleOffset(1.1);
  g1overg2erhist->GetYaxis()->SetTitleOffset(1.5);
  g1overg2erhist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g1overg2erhist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g1erhist->GetXaxis()->SetTitle("Fitted G1 uncert");
  g1erhist->GetXaxis()->SetTitleOffset(1.1);
  g1erhist->GetYaxis()->SetTitleOffset(1.5);
  g1erhist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g1erhist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  g2erhist->GetXaxis()->SetTitle("Fitted G2 uncert");
  g2erhist->GetXaxis()->SetTitleOffset(1.1);
  g2erhist->GetYaxis()->SetTitleOffset(1.5);
  g2erhist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/g2erhist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  diffFromLinG1hist->GetXaxis()->SetTitle("Normalised residual G1");
  diffFromLinG1hist->GetYaxis()->SetTitle("");
  diffFromLinG1hist->GetXaxis()->SetTitleOffset(1.1);
  diffFromLinG1hist->GetYaxis()->SetTitleOffset(1.5);
  diffFromLinG1hist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG1hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  diffFromLinG2hist->GetXaxis()->SetTitle("Normalised residual G2");
  diffFromLinG2hist->GetYaxis()->SetTitle("");
  diffFromLinG2hist->GetXaxis()->SetTitleOffset(1.1);
  diffFromLinG2hist->GetYaxis()->SetTitleOffset(1.5);
  diffFromLinG2hist->Draw();
  sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG2hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  cout << "G1 underflow " << diffFromLinG1hist->GetBinContent(0) << endl;
  cout << "G2 underflow " << diffFromLinG2hist->GetBinContent(0) << endl;
  cout << "G1 overflow " << diffFromLinG1hist->GetBinContent(101) << endl;
  cout << "G2 overflow " << diffFromLinG2hist->GetBinContent(101) << endl;

  g1overg2_isEdge->SetLineColor(kBlue);
  g1overg2_isInnerEdge->SetLineColor(kCyan);
  g1overg2_isDouble->SetLineColor(kGreen+2);
  g1overg2_isNextToDouble->SetLineColor(kRed);
  g1overg2_isQuad->SetLineColor(kOrange);

  g1overg2_isEdge->Scale(1./g1overg2_isEdge->GetEntries());
  g1overg2_isInnerEdge->Scale(1./g1overg2_isInnerEdge->GetEntries());
  g1overg2_isDouble->Scale(1./g1overg2_isDouble->GetEntries());
  g1overg2_isNextToDouble->Scale(1./g1overg2_isNextToDouble->GetEntries());
  g1overg2_isQuad->Scale(1./g1overg2_isQuad->GetEntries());
  g1overg2_isBulk->Scale(1./g1overg2_isBulk->GetEntries());

  TLegend *leg = new TLegend(0.62,0.6,0.93,0.93);
  leg->AddEntry(g1overg2_isBulk, "Normal", "l");
  leg->AddEntry(g1overg2_isDouble, "Double", "l");
  leg->AddEntry(g1overg2_isNextToDouble, "Next to D", "l");
  leg->AddEntry(g1overg2_isEdge, "Edge", "l");
  leg->AddEntry(g1overg2_isInnerEdge, "Inner E", "l");

  g1overg2_isDouble->GetXaxis()->SetTitle("G0 / G1");
  g1overg2_isDouble->GetYaxis()->SetTitle("Normalised");
  g1overg2_isDouble->GetYaxis()->SetTitleOffset(1.3);
  g1overg2_isDouble->SetMinimum(0.0);
  g1overg2_isDouble->SetMaximum(0.16);
  g1overg2_isDouble->Draw();
  g1overg2_isEdge->Draw("same");
  g1overg2_isInnerEdge->Draw("same");
  g1overg2_isNextToDouble->Draw("same");
  g1overg2_isBulk->Draw("same");
  leg->Draw("same");
  sprintf(savename,"plots/M%s/CurrentSource/g1overg2_perType.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  mapcanvas->cd();

  g1map->GetXaxis()->SetTitle("Column");
  g1map->GetYaxis()->SetTitle("Row");
  g1map->GetYaxis()->SetTitleOffset(0.7);
  if (module_str == "032") {
    g1map->GetZaxis()->SetRangeUser(-12000,-4000);
  } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008" || module_str == "044") {
    g1map->GetZaxis()->SetRangeUser(-3000,0);
  }
  g1map->Draw("colz");
  sprintf(savename,"plots/M%s/CurrentSource/g1map.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g2map->GetXaxis()->SetTitle("Column");
  g2map->GetYaxis()->SetTitle("Row");
  g2map->GetYaxis()->SetTitleOffset(0.7);
  if (module_str == "032") {
    g2map->GetZaxis()->SetRangeUser(-800,-400);
  } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008" || module_str == "044") {
    g2map->GetZaxis()->SetRangeUser(-200,-50);
  }
  g2map->Draw("colz");
  sprintf(savename,"plots/M%s/CurrentSource/g2map.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

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

  g1overg2map->GetXaxis()->SetRangeUser(0,256);
  g1overg2map->GetYaxis()->SetRangeUser(240,280);
  g1overg2map->Draw("colz");
  sprintf(savename,"plots/M%s/CurrentSource/g1overg2map_zoom.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  g1overg2ermap->GetXaxis()->SetTitle("Column");
  g1overg2ermap->GetYaxis()->SetTitle("Row");
  g1overg2ermap->GetYaxis()->SetTitleOffset(0.7);
  g1overg2ermap->GetZaxis()->SetRangeUser(0,0.5);
  g1overg2ermap->Draw("colz");
  sprintf(savename,"plots/M%s/CurrentSource/g1overg2ermap.png", 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,50);
  sprintf(savename,"plots/M%s/CurrentSource/g1ermap.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

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

  diffFromLinG1map->GetXaxis()->SetTitle("Column");
  diffFromLinG1map->GetYaxis()->SetTitle("Row");
  diffFromLinG1map->GetYaxis()->SetTitleOffset(0.7);
  diffFromLinG1map->Draw("colz");
  diffFromLinG1map->GetZaxis()->SetRangeUser(-0.04,0.04);
  sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG1map.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));
 
  diffFromLinG2map->GetXaxis()->SetTitle("Column");
  diffFromLinG2map->GetYaxis()->SetTitle("Row");
  diffFromLinG2map->GetYaxis()->SetTitleOffset(0.7);
  diffFromLinG2map->Draw("colz");
  diffFromLinG2map->GetZaxis()->SetRangeUser(-0.1,0.1);
  sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG2map.png", module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

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

}