JFCalibration/BP_DB_Comp.cpp

// to compare the results of backplane pulsing with direct beam

#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 "TFile.h"
#include "TH2F.h"
#include "TCanvas.h"

int main(int argc, char* argv[]) {
  jungfrauStyle();
  gROOT->SetBatch(1);
  gStyle->SetOptFit(11);

  if (argc != 2) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    exit(1);
  }

  string module_str = argv[1];
  char savename[128];

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

  // BP data
  sprintf(savename,"../JFModuleCalib_CuFBPCS/data/M%s/BP_ratio_M%s.root", module_str.c_str(), module_str.c_str());
  TFile* saved_file2 = new TFile((const char *)(savename),"READ");
  TH2F* g0overg1_map_BP = (TH2F*)saved_file2->Get("g0overg1_map");

  // DB data
  sprintf(savename,"data/M%s/DB_ratio_M%s.root", module_str.c_str(), module_str.c_str());
  TFile* saved_file = new TFile((const char *)(savename),"READ");
  TH2F* g0overg1_map_DB = (TH2F*)saved_file->Get("g0overg1map");
  TH2F* g0overg1er_map_DB = (TH2F*)saved_file->Get("g0overg1ermap");

  // difference and pull
  TH1F* diff_g0overg1_hist = new TH1F("diff_g0overg1_hist","",100,-20,20);
  TH2F* diff_g0overg1_map = new TH2F("diff_g0overg1_map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH1F* pull_hist = new TH1F("pull_hist","",100,-20,20);

  // scatter plots
  TH2F* scat = new TH2F("scat","",100,-50,0,100,-40,-10);
  TH2F* scat_e = new TH2F("scat_e","",100,-50,0,100,-40,-10);
  TH2F* scat_d = new TH2F("scat_d","",100,-50,0,100,-40,-10);
  TH2F* scat_ntd = new TH2F("scat_ntd","",100,-50,0,100,-40,-10);
  TH2F* scat_q = new TH2F("scat_q","",100,-50,0,100,-40,-10);
  TH2F* scat_ie = new TH2F("scat_ie","",100,-50,0,100,-40,-10);
  TH2F* scat_b = new TH2F("scat_b","",100,-50,0,100,-40,-10);

  double this_BP = 0;
  double this_DB = 0;

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

    this_BP = g0overg1_map_BP->GetBinContent((i%NC)+1,(i/NC)+1);
    this_DB = g0overg1_map_DB->GetBinContent((i%NC)+1,(i/NC)+1);

    if (this_BP != 0 && this_DB != 0) {
      diff_g0overg1_hist->Fill(this_DB - this_BP);
      diff_g0overg1_map->SetBinContent((i%NC)+1,(i/NC)+1,this_DB - this_BP);
      pull_hist->Fill((this_DB - this_BP) / g0overg1er_map_DB->GetBinContent((i%NC)+1,(i/NC)+1));

      scat->Fill(this_DB, this_BP, 1);
      if (isEdge(i)) {
	scat_e->Fill(this_DB, this_BP, 1);
      }
      if (isDouble(i)) {
	scat_d->Fill(this_DB, this_BP, 1);
      }
      if (isNextToDouble(i)) {
	scat_ntd->Fill(this_DB, this_BP, 1);
      }
      if (isQuad(i)) {
	scat_q->Fill(this_DB, this_BP, 1);
      }
      if (isInnerEdge(i)) {
	scat_ie->Fill(this_DB, this_BP, 1);
      }
      if (isBulk(i)) {
	scat_b->Fill(this_DB, this_BP, 1);
      }

    }
  }

  c1->cd();
  diff_g0overg1_hist->GetXaxis()->SetTitle("(G0 / G1)_{DB} - (G0 / G1)_{BP}");
  diff_g0overg1_hist->Draw();
  sprintf(savename,"plots/M%s/BP_DB_Comp/diff_g0overg1_hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));
  c1->SetLogy();
  diff_g0overg1_hist->Draw();
  sprintf(savename,"plots/M%s/BP_DB_Comp/diff_g0overg1_hist_log.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));
  c1->SetLogy(0);

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

  c1->cd();
  pull_hist->GetXaxis()->SetTitle("(G0 / G1)_{DB} - (G0 / G1)_{BP} / #sigma_{(G0 / G1)_{DB}}");
  pull_hist->Fit("gaus");
  pull_hist->Draw();
  sprintf(savename,"plots/M%s/BP_DB_Comp/pull_hist.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  c1->SetRightMargin(0.2);

  // all
  scat->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat->GetYaxis()->SetTitleOffset(1.3);
  scat->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // edge pixels
  scat_e->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_e->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_e->GetYaxis()->SetTitleOffset(1.3);
  scat_e->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_e.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // double pixels
  scat_d->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_d->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_d->GetYaxis()->SetTitleOffset(1.3);
  scat_d->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_d.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // next to double pixels
  scat_ntd->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_ntd->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_ntd->GetYaxis()->SetTitleOffset(1.3);
  scat_ntd->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_ntd.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // quad pixels
  scat_q->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_q->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_q->GetYaxis()->SetTitleOffset(1.3);
  scat_q->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_q.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // inner edge pixels
  scat_ie->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_ie->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_ie->GetYaxis()->SetTitleOffset(1.3);
  scat_ie->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_ie.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));

  // bulk pixels
  scat_b->GetXaxis()->SetTitle("(G0 / G1)_{DB}");
  scat_b->GetYaxis()->SetTitle("(G0 / G1)_{BP}");
  scat_b->GetYaxis()->SetTitleOffset(1.3);
  scat_b->Draw("colz");
  sprintf(savename,"plots/M%s/BP_DB_Comp/scat_b.png", module_str.c_str());
  c1->SaveAs((const char *)(savename));


}