lsbs_cert.cc

Go to the documentation of this file.

#include <Riostream.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TKey.h>
#include <TH1.h>
#include <TH2.h>
#include <TH2D.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TPaveStats.h>
#include <TText.h>
#include <TLegend.h>
#include <cstring>
#include <utility>
#include <vector>
#include <sstream>
#include <algorithm>
#include <TString.h>
#include <TColor.h>

using namespace std;

//global vars
bool debug = false;
int numlumis = -1;

int nlumis(string filename);                                            //get number of run lumisections
string runnum_str(string filename);                                     //read the run number, return in string
int getplot(string filename, string iDir, string strplot, TH1F& plot);  //read given plot
void Cleaning(vector<int>&);
string ListOut(vector<int>&);
void vector_AND(vector<int>&, vector<int>);
void lsbs_cert(string filename);

int main(int argc, char* argv[]) {
  if (argc == 2) {
    char* filename = argv[1];

    std::cout << "ready to run lsbs filename " << filename << std::endl;

    lsbs_cert(filename);

  } else {
    std::cout << "Too few arguments: " << argc << std::endl;
    return -1;
  }
  return 0;
}

void lsbs_cert(string filename) {
  void check_offset(string filename, string iDir, string plot, float limit_min, float limit_max, vector<int>&);
  void check_sigma(string filename, string iDir, string plot, float limit_err, vector<int>&);
  bool check_isgood(vector<int>&, int ls);  //check if this LS is good

  //presets
  numlumis = -1;

  float limit_x = 0.002;
  float limit_y = 0.002;
  float limit_z = 0.5;
  float limit_dx = 0.002;
  float limit_dy = 0.002;
  float limit_dz = 0.5;
  float limit_errdx = 0.002;
  float limit_errdy = 0.002;
  float limit_errdz = 0.5;

  //LS certification
  vector<int> ls_x_bad;
  vector<int> ls_y_bad;
  vector<int> ls_z_bad;

  vector<int> ls_xsc_bad;
  vector<int> ls_ysc_bad;
  vector<int> ls_zsc_bad;

  vector<int> ls_dx_bad;
  vector<int> ls_dy_bad;
  vector<int> ls_dz_bad;

  vector<int> ls_dxsc_bad;
  vector<int> ls_dysc_bad;
  vector<int> ls_dzsc_bad;

  vector<int> ls_errdx_bad;
  vector<int> ls_errdy_bad;
  vector<int> ls_errdz_bad;

  vector<int> ls_errdxsc_bad;
  vector<int> ls_errdysc_bad;
  vector<int> ls_errdzsc_bad;

  vector<int> ls_good;
  vector<int> ls_bad;

  //beamspot vs primary vertex
  check_offset(filename, "Validation", "hxLumibased PrimaryVertex-DataBase", -limit_x, limit_x, ls_x_bad);
  check_offset(filename, "Validation", "hyLumibased PrimaryVertex-DataBase", -limit_y, limit_y, ls_y_bad);
  check_offset(filename, "Validation", "hzLumibased PrimaryVertex-DataBase", -limit_z, limit_z, ls_z_bad);

  //beamspot vs scalers
  check_offset(filename, "Validation", "hxLumibased Scalers-DataBase fit", -limit_x, limit_x, ls_xsc_bad);
  check_offset(filename, "Validation", "hyLumibased Scalers-DataBase fit", -limit_y, limit_y, ls_ysc_bad);
  check_offset(filename, "Validation", "hzLumibased Scalers-DataBase fit", -limit_z, limit_z, ls_zsc_bad);

  check_offset(filename, "Debug", "hsigmaXLumibased PrimaryVertex-DataBase fit", -limit_dx, limit_dx, ls_dx_bad);
  check_offset(filename, "Debug", "hsigmaYLumibased PrimaryVertex-DataBase fit", -limit_dy, limit_dy, ls_dy_bad);
  check_offset(filename, "Debug", "hsigmaZLumibased PrimaryVertex-DataBase fit", -limit_dz, limit_dz, ls_dz_bad);

  check_offset(filename, "Validation", "hsigmaXLumibased Scalers-DataBase fit", -limit_dx, limit_dx, ls_dxsc_bad);
  check_offset(filename, "Validation", "hsigmaYLumibased Scalers-DataBase fit", -limit_dy, limit_dy, ls_dysc_bad);
  check_offset(filename, "Validation", "hsigmaZLumibased Scalers-DataBase fit", -limit_dz, limit_dz, ls_dzsc_bad);

  check_sigma(filename, "Debug", "hsigmaXLumibased PrimaryVertex-DataBase fit", limit_errdx, ls_errdx_bad);
  check_sigma(filename, "Debug", "hsigmaYLumibased PrimaryVertex-DataBase fit", limit_errdy, ls_errdy_bad);
  check_sigma(filename, "Debug", "hsigmaZLumibased PrimaryVertex-DataBase fit", limit_errdz, ls_errdz_bad);

  check_sigma(filename, "Validation", "hsigmaXLumibased Scalers-DataBase fit", limit_errdx, ls_errdxsc_bad);
  check_sigma(filename, "Validation", "hsigmaYLumibased Scalers-DataBase fit", limit_errdy, ls_errdysc_bad);
  check_sigma(filename, "Validation", "hsigmaZLumibased Scalers-DataBase fit", limit_errdz, ls_errdzsc_bad);

  //BAD LS only if bad in both histos (wrt PV, Scalers)
  vector_AND(ls_x_bad, ls_xsc_bad);
  vector_AND(ls_y_bad, ls_ysc_bad);
  vector_AND(ls_z_bad, ls_zsc_bad);
  vector_AND(ls_dx_bad, ls_dxsc_bad);
  vector_AND(ls_dy_bad, ls_dysc_bad);
  vector_AND(ls_dz_bad, ls_dzsc_bad);
  vector_AND(ls_errdx_bad, ls_errdxsc_bad);
  vector_AND(ls_errdy_bad, ls_errdysc_bad);
  vector_AND(ls_errdz_bad, ls_errdzsc_bad);

  //good LS = all LS minus BAD LS
  for (int i = 1; i <= nlumis(filename); i++) {
    if (!check_isgood(ls_x_bad, i) && !check_isgood(ls_xsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_y_bad, i) && !check_isgood(ls_ysc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_z_bad, i) && !check_isgood(ls_zsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_dx_bad, i) && !check_isgood(ls_dxsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_dy_bad, i) && !check_isgood(ls_dysc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_dz_bad, i) && !check_isgood(ls_dzsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_errdx_bad, i) && !check_isgood(ls_errdxsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_errdy_bad, i) && !check_isgood(ls_errdysc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    } else if (!check_isgood(ls_errdz_bad, i) && !check_isgood(ls_errdzsc_bad, i)) {
      ls_bad.push_back(i);
      continue;
    }

    //check also that LS is not missing!!!
    ls_good.push_back(i);
  }

  std::ofstream outfile;
  string namefile = "Certification_BS_run_" + runnum_str(filename) + ".txt";
  outfile.open(namefile.c_str());
  outfile << "Lumibased BeamSpot Calibration Certification for run " << runnum_str(filename) << ":" << endl << endl;

  char line[2000];
  sprintf(line, "    GOOD Lumisections (values within limits): %s", ListOut(ls_good).c_str());
  outfile << line << endl;

  if (!ls_bad.empty()) {
    sprintf(line, "    BAD Lumisections (values outside limits): %s", ListOut(ls_bad).c_str());
    outfile << line << endl;

    sprintf(line, "      --- histogram name ---                                --- bad lumisection list(*) ---");
    outfile << line << endl;
    sprintf(line, "      hxLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_x_bad).c_str());
    if (!ls_x_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hyLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_y_bad).c_str());
    if (!ls_y_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hzLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_z_bad).c_str());
    if (!ls_z_bad.empty())
      outfile << line << endl;

    sprintf(line, "      hsigmaXLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dx_bad).c_str());
    if (!ls_dx_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hsigmaYLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dy_bad).c_str());
    if (!ls_dy_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hsigmaZLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dz_bad).c_str());
    if (!ls_dz_bad.empty())
      outfile << line << endl;

    sprintf(line, "      hsigmaXLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdx_bad).c_str());
    if (!ls_errdx_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hsigmaYLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdy_bad).c_str());
    if (!ls_errdy_bad.empty())
      outfile << line << endl;
    sprintf(line, "      hsigmaZLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdz_bad).c_str());
    if (!ls_errdz_bad.empty())
      outfile << line << endl;

    sprintf(line, "    (*) also bad in the corresponding 'Scalers-Database fit' histograms");
    outfile << line << endl;
  }

  outfile.close();
  std::cout << "Lumibased BeamSpot Calibration Certification summary saved in " << namefile << endl;
}

void check_offset(string filename, string iDir, string plot, float limit_min, float limit_max, vector<int>& badLS) {
  TH1F checkPlot;
  if (getplot(filename, iDir, plot, checkPlot) < 0)
    return;

  //look at each LS, save the bad one
  for (int i = 1; i <= checkPlot.GetNbinsX(); i++) {
    float value = checkPlot.GetBinContent(i);
    if (value < limit_min || value > limit_max) {
      badLS.push_back((int)checkPlot.GetBinCenter(i));
    }
  }
}

void check_sigma(string filename, string iDir, string plot, float limit_err, vector<int>& badLS) {
  TH1F checkPlot;
  if (getplot(filename, iDir, plot, checkPlot) < 0)
    return;

  //look at each LS
  for (int i = 1; i <= checkPlot.GetNbinsX(); i++) {
    float value = checkPlot.GetBinError(i);
    if (value > limit_err) {
      badLS.push_back((int)checkPlot.GetBinCenter(i));
    }
  }
}

bool check_isgood(vector<int>& ls_badlist, int ls) {
  //check if this LS is found in the BAD list
  for (unsigned int i = 0; i < ls_badlist.size(); i++) {
    if (ls == ls_badlist[i])
      return false;
  }
  return true;
}

int nlumis(string filename) {
  if (numlumis > -1)
    return numlumis;

  TFile* file = TFile::Open(filename.c_str());
  if (!file->IsOpen()) {
    cerr << "Failed to open " << filename << endl;
    return -1;
  }

  string run = runnum_str(filename);
  string EventInfoDir = "DQMData/Run " + run + "/SiStrip/Run summary/EventInfo";
  TDirectory* rsEventInfoDir = dynamic_cast<TDirectory*>(file->Get(EventInfoDir.c_str()));
  rsEventInfoDir->cd();
  TIter eiKeys(rsEventInfoDir->GetListOfKeys());
  TKey* eiKey;
  while ((eiKey = dynamic_cast<TKey*>(eiKeys()))) {
    string classname(eiKey->GetClassName());
    if (classname == "TObjString") {
      string sflag = eiKey->GetName();
      string tempname = sflag.substr(sflag.find("i=") + 2);
      size_t pos1 = tempname.find('<');
      size_t pos2 = sflag.find_first_of('>');
      string detvalue = tempname.substr(0, pos1);
      string numlumisec = sflag.substr(1, pos2 - 1);
      if (numlumisec == (string) "iLumiSection") {
        numlumis = atoi(detvalue.c_str());
        break;
      }
    }
  }

  return numlumis;
}

string runnum_str(string filename) { return filename.substr(filename.find("_R000") + 5, 6); }

int getplot(string filename, string iDir, string strplot, TH1F& plot) {
  string run = runnum_str(filename);
  if (debug)
    std::cout << filename.c_str() << endl;

  TFile* file = TFile::Open(filename.c_str());
  if (!file->IsOpen()) {
    cerr << "Failed to open " << filename << endl;
    return -1;
  }

  string dir = "DQMData/Run " + run + "/AlcaBeamMonitor/Run summary/" + iDir;

  file->cd(dir.c_str());

  string theplot = strplot + ";1";
  TH1F* thisplot;
  gDirectory->GetObject(theplot.c_str(), thisplot);

  if (!thisplot) {
    std::cout << "Error: plot " << dir << "/" << theplot.c_str() << " not found!" << endl;
    return -2;
  }

  plot = *thisplot;
  thisplot = nullptr;
  delete thisplot;

  return 0;
}

void Cleaning(vector<int>& LSlist) {
  if (LSlist.empty())
    return;

  //cleaning: keep only 1st and last lumisection in the range
  int refLS = LSlist[0];
  for (unsigned int at = 1; at < LSlist.size() - 1; at++) {
    //delete LSnums in between a single continuous range
    if (refLS + 1 == LSlist[at] && LSlist[at] + 1 == LSlist[at + 1]) {
      refLS = LSlist[at];
      LSlist[at] = -1;
    } else {
      refLS = LSlist[at];
    }
  }
}

string ListOut(vector<int>& LSlist) {
  Cleaning(LSlist);

  string strout = "";
  bool rangeset = false;
  for (unsigned int at = 0; at < LSlist.size(); at++) {
    if (LSlist[at] != -1) {
      if (at > 0 && LSlist[at - 1] != -1)
        strout += ",";
      stringstream lsnum;
      lsnum << LSlist[at];
      strout += lsnum.str();
      rangeset = false;
    }
    if (LSlist[at] == -1 && !rangeset) {
      strout += "-";
      rangeset = true;
    }
  }

  return strout;
}

void vector_AND(vector<int>& bad_def, vector<int> bad_sc) {
  vector<int> temp;

  int def_size = bad_def.size();
  for (int i = 0; i < def_size; i++)
    for (unsigned int j = 0; j < bad_sc.size(); j++)
      if (bad_def[i] == bad_sc[j]) {
        temp.push_back(bad_def[i]);
        break;
      }

  bad_def = temp;
}