// to count missing and crazy calibration constants for all modules and compare

#include "../sls_detector_calibration/jungfrauCommonHeader.h"
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"

#include <fstream>
#include <dirent.h>

#include "TCanvas.h"
#include "TFile.h"
#include "TH2F.h"
#include "TGraph.h"
#include "TF1.h"

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

  jungfrauStyle();

  vector<double> modules_d;
  vector<double> hg0_missing;
  vector<double> hg0_crazy;
  vector<double> g0_missing;
  vector<double> g0_crazy;
  vector<double> g1_missing;
  vector<double> g1_crazy;
  vector<double> g2_missing;
  vector<double> g2_crazy;
  vector<double> anygain_missing_or_crazy;

  char savename[128];

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

    // check if directory exists
    sprintf(savename,"data/M%3.3d", j);    
    DIR* dir = opendir(savename);
    if (not dir) {
      continue;
    }

    int this_hg0_missing = 0;
    int this_hg0_crazy = 0;
    int this_g0_missing = 0;
    int this_g0_crazy = 0;
    int this_g1_missing = 0;
    int this_g1_crazy = 0;
    int this_g2_missing = 0;
    int this_g2_crazy = 0;
    int this_anygain_missing_or_crazy = 0;

    // CuFluo HG0 dataset
    sprintf(savename,"data/M%3.3d/CuFluo_gain_HG0_M%3.3d.root", j, j);
    TFile* FL_HG0_file = new TFile((char*)savename,"READ");
    TH2F* FL_HG0_gain_map = 0;
    if (FL_HG0_file->IsZombie()) {
      cout << "didn't find HG0 file" << endl;
    } else {
      FL_HG0_gain_map = (TH2F*)FL_HG0_file->Get("gain_ADUper1keV_2d");
    }

    // CuFluo G0 dataset
    sprintf(savename,"data/M%3.3d/CuFluo_gain_G0_M%3.3d.root", j, j);
    TFile* FL_G0_file = new TFile((char*)savename,"READ");
    TH2F* FL_G0_gain_map = 0;
    if (FL_G0_file->IsZombie()) {
      cout << "didn't find G0 file" << endl;
    } else {
      FL_G0_gain_map = (TH2F*)FL_G0_file->Get("gain_ADUper1keV_2d");
    }

    // Backplane pulsing dataset
    sprintf(savename,"data/M%3.3d/BP_ratio_M%3.3d.root", j, j);
    TFile* DB_file = new TFile((char*)savename,"READ");
    TH2F* DB_ratio_map = 0;
    if (DB_file->IsZombie()) {
      cout << "didn't find BP file" << endl;

      // look for a direct beam dataset
      sprintf(savename,"data/M%3.3d/DB_ratio_M%3.3d.root", j, j);
      DB_file = new TFile((char*)savename,"READ");

      if (DB_file->IsZombie()) {
	cout << "also didn't find DB file" << endl;

      } else {
	cout << "loading G0/G1 from DB" << endl;
	DB_ratio_map = (TH2F*)DB_file->Get("g0overg1map");
      }
    
    } else {
      cout << "loading G0/G1 from BP" << endl;
      DB_ratio_map = (TH2F*)DB_file->Get("g0overg1_map");
    }

    // Current source dataset
    sprintf(savename,"data/M%3.3d/CS_ratio_M%3.3d.root", j, j);
    TFile* CS_file = new TFile((char*)savename,"READ");
    TH2F* CS_ratio_map = 0;
    if (CS_file->IsZombie()) {
      cout << "didn't find CS file" << endl;
    } else {
      CS_ratio_map = (TH2F*)CS_file->Get("g1overg2map");
    }

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

      double this_hg0 = 0;
      double this_g0 = 0;
      double this_g0overg1 = 0;
      double this_g1 = 0;
      double this_g1overg2 = 0;
      double this_g2 = 0;
      bool this_flag = false;

      // HG0
      if (FL_HG0_gain_map) {
	this_hg0 = FL_HG0_gain_map->GetBinContent((i%NC)+1,(i/NC)+1);
	if (this_hg0 == 0) {
	  this_hg0_missing++;
	  this_flag = true;
	} else if (this_hg0 < 80 || this_hg0 > 120) {
	  this_hg0_crazy++;
	  this_flag = true;
	}
      } else {
	this_hg0_missing++;
	this_flag = true;
      }

      // G0
      if (FL_G0_gain_map) {
	this_g0 = FL_G0_gain_map->GetBinContent((i%NC)+1,(i/NC)+1);
	if (this_g0 == 0) {
	  this_g0_missing++;
	  this_flag = true;
	} else if (this_g0  < 30 || this_g0 > 48) {
	  this_g0_crazy++;
	  this_flag = true;
	}
      } else {
	this_g0_missing++;
	this_flag = true;
      }

      // G1
      if (DB_ratio_map) {
	this_g0overg1 = DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1);
	this_g1 = this_g0 / this_g0overg1;
	if (this_g0overg1 == 0) {
	  this_g1_missing++;
	  this_flag = true;
	} else if (this_g0 != 0 && (this_g1 < -2.0 || this_g1 > -0.5)) {
	  this_g1_crazy++;
	  this_flag = true;
	}
      } else {
	this_g1_missing++;
	this_flag = true;
      }

      // G2
      if (CS_ratio_map) {
	this_g1overg2 = CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1);
	this_g2 = this_g1 / this_g1overg2;
	if (this_g1overg2 == 0) {
	  this_g2_missing++;
	  this_flag = true;
	} else if (this_g0 != 0 && this_g1 != 0 && (this_g2 < -0.2 || this_g2 > -0.02)) {
	  this_g2_crazy++;
	  this_flag = true;
	}
      } else {
	this_g2_missing++;
	this_flag = true;
      }
      
      if (this_flag == true) {
	this_anygain_missing_or_crazy++;
      }

    }

    modules_d.push_back(j);
    hg0_missing.push_back(this_hg0_missing);
    hg0_crazy.push_back(this_hg0_crazy);
    g0_missing.push_back(this_g0_missing);
    g0_crazy.push_back(this_g0_crazy);
    g1_missing.push_back(this_g1_missing);
    g1_crazy.push_back(this_g1_crazy);
    g2_missing.push_back(this_g2_missing);
    g2_crazy.push_back(this_g2_crazy);
    anygain_missing_or_crazy.push_back(this_anygain_missing_or_crazy);

  }

  TCanvas *c1 = new TCanvas("c1","",150,10,800,400);
  c1->SetLeftMargin(0.12);
  c1->SetRightMargin(0.03);
  c1->SetTopMargin(0.08);
  c1->SetBottomMargin(0.15);

  c1->SetLogy();

  TF1* flat_module = new TF1("flat_module","[0]",modules_d.front(),modules_d.back());
  flat_module->SetParameter(0,524288);
  TF1* flat_chip = new TF1("flat_chip","[0]",modules_d.front(),modules_d.back());
  flat_chip->SetParameter(0,65536);
  TF1* flat_column = new TF1("flat_column","[0]",modules_d.front(),modules_d.back());
  flat_column->SetParameter(0,256);

  // HG0
  TGraph* grap_hg0_missing = new TGraph(modules_d.size(), &(modules_d[0]), &(hg0_missing[0]));
  grap_hg0_missing->SetMarkerStyle(20);
  grap_hg0_missing->GetXaxis()->SetTitle("Module");
  grap_hg0_missing->GetYaxis()->SetTitle("Pixels missing HG0 calibration");
  grap_hg0_missing->GetYaxis()->SetTitleOffset(0.9);
  grap_hg0_missing->SetMinimum(1);
  grap_hg0_missing->SetMaximum(1E6);
  grap_hg0_missing->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_hg0_missing.png");

  TGraph* grap_hg0_crazy = new TGraph(modules_d.size(), &(modules_d[0]), &(hg0_crazy[0]));
  grap_hg0_crazy->SetMarkerStyle(20);
  grap_hg0_crazy->GetXaxis()->SetTitle("Module");
  grap_hg0_crazy->GetYaxis()->SetTitle("Pixels bad HG0 calibration");
  grap_hg0_crazy->GetYaxis()->SetTitleOffset(0.9);
  grap_hg0_crazy->SetMinimum(1);
  grap_hg0_crazy->SetMaximum(1E6);
  grap_hg0_crazy->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_hg0_crazy.png");

  // G0
  TGraph* grap_g0_missing = new TGraph(modules_d.size(), &(modules_d[0]), &(g0_missing[0]));
  grap_g0_missing->SetMarkerStyle(20);
  grap_g0_missing->GetXaxis()->SetTitle("Module");
  grap_g0_missing->GetYaxis()->SetTitle("Pixels missing G0 calibration");
  grap_g0_missing->GetYaxis()->SetTitleOffset(0.9);
  grap_g0_missing->SetMinimum(1);
  grap_g0_missing->SetMaximum(1E6);
  grap_g0_missing->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g0_missing.png");

  TGraph* grap_g0_crazy = new TGraph(modules_d.size(), &(modules_d[0]), &(g0_crazy[0]));
  grap_g0_crazy->SetMarkerStyle(20);
  grap_g0_crazy->GetXaxis()->SetTitle("Module");
  grap_g0_crazy->GetYaxis()->SetTitle("Pixels bad G0 calibration");
  grap_g0_crazy->GetYaxis()->SetTitleOffset(0.9);
  grap_g0_crazy->SetMinimum(1);
  grap_g0_crazy->SetMaximum(1E6);
  grap_g0_crazy->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g0_crazy.png");

  // G1
  TGraph* grap_g1_missing = new TGraph(modules_d.size(), &(modules_d[0]), &(g1_missing[0]));
  grap_g1_missing->SetMarkerStyle(20);
  grap_g1_missing->GetXaxis()->SetTitle("Module");
  grap_g1_missing->GetYaxis()->SetTitle("Pixels missing G1 calibration");
  grap_g1_missing->GetYaxis()->SetTitleOffset(0.9);
  grap_g1_missing->SetMinimum(1);
  grap_g1_missing->SetMaximum(1E6);
  grap_g1_missing->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g1_missing.png");

  TGraph* grap_g1_crazy = new TGraph(modules_d.size(), &(modules_d[0]), &(g1_crazy[0]));
  grap_g1_crazy->SetMarkerStyle(20);
  grap_g1_crazy->GetXaxis()->SetTitle("Module");
  grap_g1_crazy->GetYaxis()->SetTitle("Pixels bad G1 calibration");
  grap_g1_crazy->GetYaxis()->SetTitleOffset(0.9);
  grap_g1_crazy->SetMinimum(1);
  grap_g1_crazy->SetMaximum(1E6);
  grap_g1_crazy->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g1_crazy.png");

  // G2
  TGraph* grap_g2_missing = new TGraph(modules_d.size(), &(modules_d[0]), &(g2_missing[0]));
  grap_g2_missing->SetMarkerStyle(20);
  grap_g2_missing->GetXaxis()->SetTitle("Module");
  grap_g2_missing->GetYaxis()->SetTitle("Pixels missing G2 calibration");
  grap_g2_missing->GetYaxis()->SetTitleOffset(0.9);
  grap_g2_missing->SetMinimum(1);
  grap_g2_missing->SetMaximum(1E6);
  grap_g2_missing->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g2_missing.png");

  TGraph* grap_g2_crazy = new TGraph(modules_d.size(), &(modules_d[0]), &(g2_crazy[0]));
  grap_g2_crazy->SetMarkerStyle(20);
  grap_g2_crazy->GetXaxis()->SetTitle("Module");
  grap_g2_crazy->GetYaxis()->SetTitle("Pixels bad G2 calibration");
  grap_g2_crazy->GetYaxis()->SetTitleOffset(0.9);
  grap_g2_crazy->SetMinimum(1);
  grap_g2_crazy->SetMaximum(1E6);
  grap_g2_crazy->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_g2_crazy.png");

  // any
  TGraph* grap_any = new TGraph(modules_d.size(), &(modules_d[0]), &(anygain_missing_or_crazy[0]));
  grap_any->SetMarkerStyle(20);
  grap_any->GetXaxis()->SetTitle("Module");
  grap_any->GetYaxis()->SetTitle("Pixels anything wrong");
  grap_any->GetYaxis()->SetTitleOffset(0.9);
  grap_any->SetMinimum(1);
  grap_any->SetMaximum(1E6);
  grap_any->Draw("AP");
  flat_module->Draw("same");
  flat_chip->Draw("same");
  flat_column->Draw("same");
  c1->SaveAs("plots/BadPixels_anythingwrong.png");

}