SiStripMonitorClient/bin/lsbs_cert.cc

0001 #include <Riostream.h>
0002 #include <TDirectory.h>
0003 #include <TFile.h>
0004 #include <TROOT.h>
0005 #include <TStyle.h>
0006 #include <TKey.h>
0007 #include <TH1.h>
0008 #include <TH2.h>
0009 #include <TH2D.h>
0010 #include <TCanvas.h>
0011 #include <TGraph.h>
0012 #include <TPaveStats.h>
0013 #include <TText.h>
0014 #include <TLegend.h>
0015 #include <cstring>
0016 #include <utility>
0017 #include <vector>
0018 #include <sstream>
0019 #include <algorithm>
0020 #include <TString.h>
0021 #include <TColor.h>
0022
0023 using namespace std;
0024
0025 //global vars
0026 bool debug = false;
0027 int numlumis = -1;
0028
0029 int nlumis(string filename);                                            //get number of run lumisections
0030 string runnum_str(string filename);                                     //read the run number, return in string
0031 int getplot(string filename, string iDir, string strplot, TH1F& plot);  //read given plot
0032 void Cleaning(vector<int>&);
0033 string ListOut(vector<int>&);
0034 void vector_AND(vector<int>&, vector<int>);
0035 void lsbs_cert(string filename);
0036
0037 int main(int argc, char* argv[]) {
0038   if (argc == 2) {
0039     char* filename = argv[1];
0040
0041     std::cout << "ready to run lsbs filename " << filename << std::endl;
0042
0043     lsbs_cert(filename);
0044
0045   } else {
0046     std::cout << "Too few arguments: " << argc << std::endl;
0047     return -1;
0048   }
0049   return 0;
0050 }
0051
0052 void lsbs_cert(string filename) {
0053   void check_offset(string filename, string iDir, string plot, float limit_min, float limit_max, vector<int>&);
0054   void check_sigma(string filename, string iDir, string plot, float limit_err, vector<int>&);
0055   bool check_isgood(vector<int>&, int ls);  //check if this LS is good
0056
0057   //presets
0058   numlumis = -1;
0059
0060   float limit_x = 0.002;
0061   float limit_y = 0.002;
0062   float limit_z = 0.5;
0063   float limit_dx = 0.002;
0064   float limit_dy = 0.002;
0065   float limit_dz = 0.5;
0066   float limit_errdx = 0.002;
0067   float limit_errdy = 0.002;
0068   float limit_errdz = 0.5;
0069
0070   //LS certification
0071   vector<int> ls_x_bad;
0072   vector<int> ls_y_bad;
0073   vector<int> ls_z_bad;
0074
0075   vector<int> ls_xsc_bad;
0076   vector<int> ls_ysc_bad;
0077   vector<int> ls_zsc_bad;
0078
0079   vector<int> ls_dx_bad;
0080   vector<int> ls_dy_bad;
0081   vector<int> ls_dz_bad;
0082
0083   vector<int> ls_dxsc_bad;
0084   vector<int> ls_dysc_bad;
0085   vector<int> ls_dzsc_bad;
0086
0087   vector<int> ls_errdx_bad;
0088   vector<int> ls_errdy_bad;
0089   vector<int> ls_errdz_bad;
0090
0091   vector<int> ls_errdxsc_bad;
0092   vector<int> ls_errdysc_bad;
0093   vector<int> ls_errdzsc_bad;
0094
0095   vector<int> ls_good;
0096   vector<int> ls_bad;
0097
0098   //beamspot vs primary vertex
0099   check_offset(filename, "Validation", "hxLumibased PrimaryVertex-DataBase", -limit_x, limit_x, ls_x_bad);
0100   check_offset(filename, "Validation", "hyLumibased PrimaryVertex-DataBase", -limit_y, limit_y, ls_y_bad);
0101   check_offset(filename, "Validation", "hzLumibased PrimaryVertex-DataBase", -limit_z, limit_z, ls_z_bad);
0102
0103   //beamspot vs scalers
0104   check_offset(filename, "Validation", "hxLumibased Scalers-DataBase fit", -limit_x, limit_x, ls_xsc_bad);
0105   check_offset(filename, "Validation", "hyLumibased Scalers-DataBase fit", -limit_y, limit_y, ls_ysc_bad);
0106   check_offset(filename, "Validation", "hzLumibased Scalers-DataBase fit", -limit_z, limit_z, ls_zsc_bad);
0107
0108   check_offset(filename, "Debug", "hsigmaXLumibased PrimaryVertex-DataBase fit", -limit_dx, limit_dx, ls_dx_bad);
0109   check_offset(filename, "Debug", "hsigmaYLumibased PrimaryVertex-DataBase fit", -limit_dy, limit_dy, ls_dy_bad);
0110   check_offset(filename, "Debug", "hsigmaZLumibased PrimaryVertex-DataBase fit", -limit_dz, limit_dz, ls_dz_bad);
0111
0112   check_offset(filename, "Validation", "hsigmaXLumibased Scalers-DataBase fit", -limit_dx, limit_dx, ls_dxsc_bad);
0113   check_offset(filename, "Validation", "hsigmaYLumibased Scalers-DataBase fit", -limit_dy, limit_dy, ls_dysc_bad);
0114   check_offset(filename, "Validation", "hsigmaZLumibased Scalers-DataBase fit", -limit_dz, limit_dz, ls_dzsc_bad);
0115
0116   check_sigma(filename, "Debug", "hsigmaXLumibased PrimaryVertex-DataBase fit", limit_errdx, ls_errdx_bad);
0117   check_sigma(filename, "Debug", "hsigmaYLumibased PrimaryVertex-DataBase fit", limit_errdy, ls_errdy_bad);
0118   check_sigma(filename, "Debug", "hsigmaZLumibased PrimaryVertex-DataBase fit", limit_errdz, ls_errdz_bad);
0119
0120   check_sigma(filename, "Validation", "hsigmaXLumibased Scalers-DataBase fit", limit_errdx, ls_errdxsc_bad);
0121   check_sigma(filename, "Validation", "hsigmaYLumibased Scalers-DataBase fit", limit_errdy, ls_errdysc_bad);
0122   check_sigma(filename, "Validation", "hsigmaZLumibased Scalers-DataBase fit", limit_errdz, ls_errdzsc_bad);
0123
0124   //BAD LS only if bad in both histos (wrt PV, Scalers)
0125   vector_AND(ls_x_bad, ls_xsc_bad);
0126   vector_AND(ls_y_bad, ls_ysc_bad);
0127   vector_AND(ls_z_bad, ls_zsc_bad);
0128   vector_AND(ls_dx_bad, ls_dxsc_bad);
0129   vector_AND(ls_dy_bad, ls_dysc_bad);
0130   vector_AND(ls_dz_bad, ls_dzsc_bad);
0131   vector_AND(ls_errdx_bad, ls_errdxsc_bad);
0132   vector_AND(ls_errdy_bad, ls_errdysc_bad);
0133   vector_AND(ls_errdz_bad, ls_errdzsc_bad);
0134
0135   //good LS = all LS minus BAD LS
0136   for (int i = 1; i <= nlumis(filename); i++) {
0137     if (!check_isgood(ls_x_bad, i) && !check_isgood(ls_xsc_bad, i)) {
0138       ls_bad.push_back(i);
0139       continue;
0140     } else if (!check_isgood(ls_y_bad, i) && !check_isgood(ls_ysc_bad, i)) {
0141       ls_bad.push_back(i);
0142       continue;
0143     } else if (!check_isgood(ls_z_bad, i) && !check_isgood(ls_zsc_bad, i)) {
0144       ls_bad.push_back(i);
0145       continue;
0146     } else if (!check_isgood(ls_dx_bad, i) && !check_isgood(ls_dxsc_bad, i)) {
0147       ls_bad.push_back(i);
0148       continue;
0149     } else if (!check_isgood(ls_dy_bad, i) && !check_isgood(ls_dysc_bad, i)) {
0150       ls_bad.push_back(i);
0151       continue;
0152     } else if (!check_isgood(ls_dz_bad, i) && !check_isgood(ls_dzsc_bad, i)) {
0153       ls_bad.push_back(i);
0154       continue;
0155     } else if (!check_isgood(ls_errdx_bad, i) && !check_isgood(ls_errdxsc_bad, i)) {
0156       ls_bad.push_back(i);
0157       continue;
0158     } else if (!check_isgood(ls_errdy_bad, i) && !check_isgood(ls_errdysc_bad, i)) {
0159       ls_bad.push_back(i);
0160       continue;
0161     } else if (!check_isgood(ls_errdz_bad, i) && !check_isgood(ls_errdzsc_bad, i)) {
0162       ls_bad.push_back(i);
0163       continue;
0164     }
0165
0166     //check also that LS is not missing!!!
0167     ls_good.push_back(i);
0168   }
0169
0170   std::ofstream outfile;
0171   string namefile = "Certification_BS_run_" + runnum_str(filename) + ".txt";
0172   outfile.open(namefile.c_str());
0173   outfile << "Lumibased BeamSpot Calibration Certification for run " << runnum_str(filename) << ":" << endl << endl;
0174
0175   char line[2000];
0176   sprintf(line, "    GOOD Lumisections (values within limits): %s", ListOut(ls_good).c_str());
0177   outfile << line << endl;
0178
0179   if (!ls_bad.empty()) {
0180     sprintf(line, "    BAD Lumisections (values outside limits): %s", ListOut(ls_bad).c_str());
0181     outfile << line << endl;
0182
0183     sprintf(line, "      --- histogram name ---                                --- bad lumisection list(*) ---");
0184     outfile << line << endl;
0185     sprintf(line, "      hxLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_x_bad).c_str());
0186     if (!ls_x_bad.empty())
0187       outfile << line << endl;
0188     sprintf(line, "      hyLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_y_bad).c_str());
0189     if (!ls_y_bad.empty())
0190       outfile << line << endl;
0191     sprintf(line, "      hzLumibased PrimaryVertex-DataBase (mean):            %s", ListOut(ls_z_bad).c_str());
0192     if (!ls_z_bad.empty())
0193       outfile << line << endl;
0194
0195     sprintf(line, "      hsigmaXLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dx_bad).c_str());
0196     if (!ls_dx_bad.empty())
0197       outfile << line << endl;
0198     sprintf(line, "      hsigmaYLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dy_bad).c_str());
0199     if (!ls_dy_bad.empty())
0200       outfile << line << endl;
0201     sprintf(line, "      hsigmaZLumibased PrimaryVertex-DataBase fit (mean):   %s", ListOut(ls_dz_bad).c_str());
0202     if (!ls_dz_bad.empty())
0203       outfile << line << endl;
0204
0205     sprintf(line, "      hsigmaXLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdx_bad).c_str());
0206     if (!ls_errdx_bad.empty())
0207       outfile << line << endl;
0208     sprintf(line, "      hsigmaYLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdy_bad).c_str());
0209     if (!ls_errdy_bad.empty())
0210       outfile << line << endl;
0211     sprintf(line, "      hsigmaZLumibased PrimaryVertex-DataBase fit (error):  %s", ListOut(ls_errdz_bad).c_str());
0212     if (!ls_errdz_bad.empty())
0213       outfile << line << endl;
0214
0215     sprintf(line, "    (*) also bad in the corresponding 'Scalers-Database fit' histograms");
0216     outfile << line << endl;
0217   }
0218
0219   outfile.close();
0220   std::cout << "Lumibased BeamSpot Calibration Certification summary saved in " << namefile << endl;
0221 }
0222
0223 void check_offset(string filename, string iDir, string plot, float limit_min, float limit_max, vector<int>& badLS) {
0224   TH1F checkPlot;
0225   if (getplot(filename, iDir, plot, checkPlot) < 0)
0226     return;
0227
0228   //look at each LS, save the bad one
0229   for (int i = 1; i <= checkPlot.GetNbinsX(); i++) {
0230     float value = checkPlot.GetBinContent(i);
0231     if (value < limit_min || value > limit_max) {
0232       badLS.push_back((int)checkPlot.GetBinCenter(i));
0233     }
0234   }
0235 }
0236
0237 void check_sigma(string filename, string iDir, string plot, float limit_err, vector<int>& badLS) {
0238   TH1F checkPlot;
0239   if (getplot(filename, iDir, plot, checkPlot) < 0)
0240     return;
0241
0242   //look at each LS
0243   for (int i = 1; i <= checkPlot.GetNbinsX(); i++) {
0244     float value = checkPlot.GetBinError(i);
0245     if (value > limit_err) {
0246       badLS.push_back((int)checkPlot.GetBinCenter(i));
0247     }
0248   }
0249 }
0250
0251 bool check_isgood(vector<int>& ls_badlist, int ls) {
0252   //check if this LS is found in the BAD list
0253   for (unsigned int i = 0; i < ls_badlist.size(); i++) {
0254     if (ls == ls_badlist[i])
0255       return false;
0256   }
0257   return true;
0258 }
0259
0260 int nlumis(string filename) {
0261   if (numlumis > -1)
0262     return numlumis;
0263
0264   TFile* file = TFile::Open(filename.c_str());
0265   if (!file->IsOpen()) {
0266     cerr << "Failed to open " << filename << endl;
0267     return -1;
0268   }
0269
0270   string run = runnum_str(filename);
0271   string EventInfoDir = "DQMData/Run " + run + "/SiStrip/Run summary/EventInfo";
0272   TDirectory* rsEventInfoDir = dynamic_cast<TDirectory*>(file->Get(EventInfoDir.c_str()));
0273   rsEventInfoDir->cd();
0274   TIter eiKeys(rsEventInfoDir->GetListOfKeys());
0275   TKey* eiKey;
0276   while ((eiKey = dynamic_cast<TKey*>(eiKeys()))) {
0277     string classname(eiKey->GetClassName());
0278     if (classname == "TObjString") {
0279       string sflag = eiKey->GetName();
0280       string tempname = sflag.substr(sflag.find("i=") + 2);
0281       size_t pos1 = tempname.find('<');
0282       size_t pos2 = sflag.find_first_of('>');
0283       string detvalue = tempname.substr(0, pos1);
0284       string numlumisec = sflag.substr(1, pos2 - 1);
0285       if (numlumisec == (string) "iLumiSection") {
0286         numlumis = atoi(detvalue.c_str());
0287         break;
0288       }
0289     }
0290   }
0291
0292   return numlumis;
0293 }
0294
0295 string runnum_str(string filename) { return filename.substr(filename.find("_R000") + 5, 6); }
0296
0297 int getplot(string filename, string iDir, string strplot, TH1F& plot) {
0298   string run = runnum_str(filename);
0299   if (debug)
0300     std::cout << filename.c_str() << endl;
0301
0302   TFile* file = TFile::Open(filename.c_str());
0303   if (!file->IsOpen()) {
0304     cerr << "Failed to open " << filename << endl;
0305     return -1;
0306   }
0307
0308   string dir = "DQMData/Run " + run + "/AlcaBeamMonitor/Run summary/" + iDir;
0309
0310   file->cd(dir.c_str());
0311
0312   string theplot = strplot + ";1";
0313   TH1F* thisplot;
0314   gDirectory->GetObject(theplot.c_str(), thisplot);
0315
0316   if (!thisplot) {
0317     std::cout << "Error: plot " << dir << "/" << theplot.c_str() << " not found!" << endl;
0318     return -2;
0319   }
0320
0321   plot = *thisplot;
0322   thisplot = nullptr;
0323   delete thisplot;
0324
0325   return 0;
0326 }
0327
0328 void Cleaning(vector<int>& LSlist) {
0329   if (LSlist.empty())
0330     return;
0331
0332   //cleaning: keep only 1st and last lumisection in the range
0333   int refLS = LSlist[0];
0334   for (unsigned int at = 1; at < LSlist.size() - 1; at++) {
0335     //delete LSnums in between a single continuous range
0336     if (refLS + 1 == LSlist[at] && LSlist[at] + 1 == LSlist[at + 1]) {
0337       refLS = LSlist[at];
0338       LSlist[at] = -1;
0339     } else {
0340       refLS = LSlist[at];
0341     }
0342   }
0343 }
0344
0345 string ListOut(vector<int>& LSlist) {
0346   Cleaning(LSlist);
0347
0348   string strout = "";
0349   bool rangeset = false;
0350   for (unsigned int at = 0; at < LSlist.size(); at++) {
0351     if (LSlist[at] != -1) {
0352       if (at > 0 && LSlist[at - 1] != -1)
0353         strout += ",";
0354       stringstream lsnum;
0355       lsnum << LSlist[at];
0356       strout += lsnum.str();
0357       rangeset = false;
0358     }
0359     if (LSlist[at] == -1 && !rangeset) {
0360       strout += "-";
0361       rangeset = true;
0362     }
0363   }
0364
0365   return strout;
0366 }
0367
0368 void vector_AND(vector<int>& bad_def, vector<int> bad_sc) {
0369   vector<int> temp;
0370
0371   int def_size = bad_def.size();
0372   for (int i = 0; i < def_size; i++)
0373     for (unsigned int j = 0; j < bad_sc.size(); j++)
0374       if (bad_def[i] == bad_sc[j]) {
0375         temp.push_back(bad_def[i]);
0376         break;
0377       }
0378
0379   bad_def = temp;
0380 }