JFCalibration/JFMC_ModuleComp.cpp

// 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 <stdexcept>

#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];
};

std::string stdoutOfCommand(const char* cmd) {
    char buffer[128];
    std::string result = "";
    FILE* pipe = popen(cmd, "r");
    if (!pipe) throw std::runtime_error("popen() failed!");
    try {
        while (fgets(buffer, sizeof buffer, pipe) != NULL) {
            result += buffer;
        }
    } catch (...) {
        pclose(pipe);
        throw;
    }
    pclose(pipe);
    return result;
}

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

  jungfrauStyle();

  GainMaps *mapsObject = new GainMaps;
  
  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 (int j = 0; j < 300; 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
    sprintf(savename,"ls -1 data/M%3.3d/gainMaps_M%3.3d* | tail -1", j, j);
    string mostRecentGainFile = stdoutOfCommand(savename);
    mostRecentGainFile.erase(std::remove(mostRecentGainFile.begin(), mostRecentGainFile.end(), '\n'), mostRecentGainFile.end());

    fstream infile;
    infile.open(mostRecentGainFile.c_str(), ios::in | ios::binary);
    if (not infile.is_open()) {
      continue;
    }
    cout << "found " << mostRecentGainFile << endl;

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