// to extract the calibration constants of all modules from the gain maps and compare

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

#include <fstream> 

#include "TH1F.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraphErrors.h"

struct GainMaps
{
  double g0VALs[NCH];
  double g1VALs[NCH];
  double g2VALs[NCH];
  double hg0VALs[NCH];
};


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

  jungfrauStyle();

  GainMaps *mapsObject = new GainMaps;

  vector<int> modules;
  modules.push_back(6);
  modules.push_back(8);
  modules.push_back(21);
  modules.push_back(22);
  modules.push_back(27);
  modules.push_back(28);
  modules.push_back(32);
  modules.push_back(35);
  modules.push_back(39);
  modules.push_back(40);
  modules.push_back(44);
  modules.push_back(45);
  modules.push_back(48);
  modules.push_back(50);
  modules.push_back(55);
  modules.push_back(57);
  modules.push_back(59);
  modules.push_back(62);
  modules.push_back(64);
  modules.push_back(65);
  modules.push_back(66);
  modules.push_back(67);
  modules.push_back(68);
  modules.push_back(69);
  modules.push_back(72);
  modules.push_back(74);
  modules.push_back(76);
  modules.push_back(79);
  modules.push_back(80);
  modules.push_back(86);
  modules.push_back(88);
  modules.push_back(94);
  modules.push_back(100);
  modules.push_back(102);
  modules.push_back(103);
  modules.push_back(104);
  modules.push_back(108);
  modules.push_back(109);
  modules.push_back(110);
  modules.push_back(113);
  modules.push_back(119);
  modules.push_back(121);
  modules.push_back(124);
  modules.push_back(126);
  modules.push_back(131);
  modules.push_back(204);
  modules.push_back(205);
  modules.push_back(208);
  modules.push_back(212);
  modules.push_back(222);
  modules.push_back(226);
  
  vector<double> modules_d;
  vector<double> module_errs;
  vector<double> g0_means;
  vector<double> g0_rmss;
  vector<double> g1_means;
  vector<double> g1_rmss;
  vector<double> g2_means;
  vector<double> g2_rmss;
  vector<double> hg0_means;
  vector<double> hg0_rmss;

  char savename[128];

  TH1F* g0hist = new TH1F("g0hist","",100,30,50);
  TH1F* g1hist = new TH1F("g1hist","",100,-3,0);
  TH1F* g2hist = new TH1F("g2hist","",100,-0.2,0);
  TH1F* hg0hist = new TH1F("hg0hist","",100,80,120);

  for (size_t j = 0; j < modules.size(); j++) {

    // clear gain maps
    for (int i = 0; i < NCH; i++) {
      mapsObject->g0VALs[i] = 0;
      mapsObject->g1VALs[i] = 0;
      mapsObject->g2VALs[i] = 0;
      mapsObject->hg0VALs[i] = 0;
    }

    // get gain maps
    fstream infile;
    sprintf(savename,"data/M%3.3d/gainMaps_M%3.3d.bin", modules[j], modules[j]);
    infile.open((char*)savename, ios::in | ios::binary);
    
    infile.read((char*)mapsObject->g0VALs, sizeof(mapsObject->g0VALs));
    infile.read((char*)mapsObject->g1VALs, sizeof(mapsObject->g1VALs));
    infile.read((char*)mapsObject->g2VALs, sizeof(mapsObject->g2VALs));
    infile.read((char*)mapsObject->hg0VALs, sizeof(mapsObject->hg0VALs));

    infile.close();

    for (int i = 0; i < NCH; i++) {
      g0hist->Fill(mapsObject->g0VALs[i]);
      g1hist->Fill(mapsObject->g1VALs[i]);
      g2hist->Fill(mapsObject->g2VALs[i]);
      hg0hist->Fill(mapsObject->hg0VALs[i]);
    }

    modules_d.push_back(modules[j]);
    module_errs.push_back(0.);
    g0_means.push_back(g0hist->GetMean());
    g0_rmss.push_back(g0hist->GetRMS());
    g1_means.push_back(g1hist->GetMean());
    g1_rmss.push_back(g1hist->GetRMS());
    g2_means.push_back(g2hist->GetMean());
    g2_rmss.push_back(g2hist->GetRMS());
    hg0_means.push_back(hg0hist->GetMean());
    hg0_rmss.push_back(hg0hist->GetRMS());

    g0hist->Reset();
    g1hist->Reset();
    g2hist->Reset();
    hg0hist->Reset();
    
  }

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

  TGraphErrors* grap_g0 = new TGraphErrors(modules_d.size(), &(modules_d[0]), &(g0_means[0]), &(module_errs[0]), &(g0_rmss[0]));
  grap_g0->SetMarkerStyle(20);
  grap_g0->GetXaxis()->SetTitle("Module");
  grap_g0->GetYaxis()->SetTitle("Mean G0 [ADU/keV]");
  grap_g0->GetYaxis()->SetTitleOffset(0.9);
  grap_g0->SetMinimum(38);
  grap_g0->SetMaximum(44);
  grap_g0->Draw("AP");
  c1->SaveAs("plots/ModuleSummary_g0.png");

  TGraphErrors* grap_g1 = new TGraphErrors(modules_d.size(), &(modules_d[0]), &(g1_means[0]), &(module_errs[0]), &(g1_rmss[0]));
  grap_g1->SetMarkerStyle(20);
  grap_g1->GetXaxis()->SetTitle("Module");
  grap_g1->GetYaxis()->SetTitle("Mean G1 [ADU/keV]");
  grap_g1->GetYaxis()->SetTitleOffset(0.9);
  grap_g1->SetMinimum(-1.6);
  grap_g1->SetMaximum(-1.2);
  grap_g1->Draw("AP");
  c1->SaveAs("plots/ModuleSummary_g1.png");

  TGraphErrors* grap_g2 = new TGraphErrors(modules_d.size(), &(modules_d[0]), &(g2_means[0]), &(module_errs[0]), &(g2_rmss[0]));
  grap_g2->SetMarkerStyle(20);
  grap_g2->GetXaxis()->SetTitle("Module");
  grap_g2->GetYaxis()->SetTitle("Mean G2 [ADU/keV]");
  grap_g2->GetYaxis()->SetTitleOffset(0.9);
  grap_g2->SetMinimum(-0.13);
  grap_g2->SetMaximum(-0.09);
  grap_g2->Draw("AP");
  c1->SaveAs("plots/ModuleSummary_g2.png");

  TGraphErrors* grap_hg0 = new TGraphErrors(modules_d.size(), &(modules_d[0]), &(hg0_means[0]), &(module_errs[0]), &(hg0_rmss[0]));
  grap_hg0->SetMarkerStyle(20);
  grap_hg0->GetXaxis()->SetTitle("Module");
  grap_hg0->GetYaxis()->SetTitle("Mean HG0 [ADU/keV]");
  grap_hg0->GetYaxis()->SetTitleOffset(0.9);
  grap_hg0->SetMinimum(94);
  grap_hg0->SetMaximum(106);
  grap_hg0->Draw("AP");
  c1->SaveAs("plots/ModuleSummary_hg0.png");
  
}