Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:08:39

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 int numlumis = -1;
0027 
0028 int nlumis(string filename);         //get number of run lumisections
0029 string runnum_str(string filename);  //read the run number, return in string
0030 void ls_cert(float threshold_pixel, float threshold, string filename);
0031 
0032 int main(int argc, char* argv[]) {
0033   if (argc == 4) {
0034     char* cpixel = argv[1];
0035     char* cthr = argv[2];
0036     char* filename = argv[3];
0037 
0038     float threshold_pixel = 0;
0039     sscanf(cpixel, "%f", &threshold_pixel);
0040     float threshold = 0;
0041     sscanf(cthr, "%f", &threshold);
0042 
0043     std::cout << "ready to run ls_cert: pixel thr  " << threshold_pixel << " threshold " << threshold << " filename "
0044               << filename << std::endl;
0045 
0046     ls_cert(threshold_pixel, threshold, filename);
0047 
0048   } else {
0049     std::cout << "Too few arguments: " << argc << std::endl;
0050     return -1;
0051   }
0052   return 0;
0053 }
0054 
0055 void ls_cert(float threshold_pixel, float threshold, string filename) {
0056   void ls_cert_type(string iDir,
0057                     float threshold,
0058                     string filename,
0059                     vector<string>&,
0060                     vector<pair<string, vector<float> > >&,
0061                     vector<pair<string, vector<float> > >&,
0062                     vector<pair<string, vector<float> > >&);
0063   void cert_plot(float threshold_pixel,
0064                  float threshold,
0065                  string filename,
0066                  vector<string>&,
0067                  vector<string>&,
0068                  vector<string>&,
0069                  vector<pair<string, vector<float> > >&,
0070                  vector<pair<string, vector<float> > >&,
0071                  vector<pair<string, vector<float> > >&);
0072 
0073   //presets
0074   numlumis = -1;
0075 
0076   //certifications
0077   vector<string> cert_strip;
0078   vector<string> cert_track;
0079   vector<string> cert_pixel;
0080 
0081   //good lumisections
0082   vector<pair<string, vector<float> > > gLS_strip;
0083   vector<pair<string, vector<float> > > gLS_track;
0084   vector<pair<string, vector<float> > > gLS_pixel;
0085 
0086   //bad lumisections
0087   vector<pair<string, vector<float> > > bLS_strip;
0088   vector<pair<string, vector<float> > > bLS_track;
0089   vector<pair<string, vector<float> > > bLS_pixel;
0090 
0091   //missing lumisections
0092   vector<pair<string, vector<float> > > mLS_strip;
0093   vector<pair<string, vector<float> > > mLS_track;
0094   vector<pair<string, vector<float> > > mLS_pixel;
0095 
0096   ls_cert_type("SiStrip", threshold, filename, cert_strip, gLS_strip, bLS_strip, mLS_strip);
0097   ls_cert_type("Tracking", threshold, filename, cert_track, gLS_track, bLS_track, mLS_track);
0098   ls_cert_type("Pixel", threshold_pixel, filename, cert_pixel, gLS_pixel, bLS_pixel, mLS_pixel);
0099 
0100   std::ofstream outfile;
0101   string namefile = "Certification_run_" + runnum_str(filename) + ".txt";
0102   outfile.open(namefile.c_str());
0103   outfile << "Lumisection Certification (GOOD: >= " << threshold_pixel << " [Pixel]; >= " << threshold
0104           << " [SiStrip,Tracking] "
0105           << ", otherwise BAD):" << endl
0106           << endl;
0107   outfile << "GOOD Lumisections:" << endl;
0108   char line[200];
0109   for (int ityp = 0; ityp < 4; ityp++) {
0110     sprintf(line, " Pixel    %*sSummary: %s", 13, cert_pixel[ityp].c_str(), gLS_pixel[ityp].first.c_str());
0111     outfile << line << endl;
0112   }
0113   for (int ityp = 0; ityp < 4; ityp++) {
0114     sprintf(line, " SiStrip  %*sSummary: %s", 13, cert_strip[ityp].c_str(), gLS_strip[ityp].first.c_str());
0115     outfile << line << endl;
0116   }
0117   for (int ityp = 0; ityp < 1; ityp++) {
0118     sprintf(line, " Tracking %*sSummary: %s", 13, cert_track[ityp].c_str(), gLS_track[ityp].first.c_str());
0119     outfile << line << endl;
0120   }
0121 
0122   outfile << "\nBAD Lumisections:" << endl;
0123   for (int ityp = 0; ityp < 4; ityp++) {
0124     sprintf(line, " Pixel    %*sSummary: %s", 13, cert_pixel[ityp].c_str(), bLS_pixel[ityp].first.c_str());
0125     outfile << line << endl;
0126   }
0127   for (int ityp = 0; ityp < 4; ityp++) {
0128     sprintf(line, " SiStrip  %*sSummary: %s", 13, cert_strip[ityp].c_str(), bLS_strip[ityp].first.c_str());
0129     outfile << line << endl;
0130   }
0131   for (int ityp = 0; ityp < 1; ityp++) {
0132     sprintf(line, " Tracking %*sSummary: %s", 13, cert_track[ityp].c_str(), bLS_track[ityp].first.c_str());
0133     outfile << line << endl;
0134   }
0135 
0136   outfile << "\nMISSING Lumisections:" << endl;
0137   for (int ityp = 0; ityp < 4; ityp++) {
0138     sprintf(line, " Pixel    %*sSummary: %s", 13, cert_pixel[ityp].c_str(), mLS_pixel[ityp].first.c_str());
0139     outfile << line << endl;
0140   }
0141   for (int ityp = 0; ityp < 4; ityp++) {
0142     sprintf(line, " SiStrip  %*sSummary: %s", 13, cert_strip[ityp].c_str(), mLS_strip[ityp].first.c_str());
0143     outfile << line << endl;
0144   }
0145   for (int ityp = 0; ityp < 1; ityp++) {
0146     sprintf(line, " Tracking %*sSummary: %s", 13, cert_track[ityp].c_str(), mLS_track[ityp].first.c_str());
0147     outfile << line << endl;
0148   }
0149 
0150   outfile.close();
0151   std::cout << "Lumisection Certification summary saved in " << namefile << std::endl;
0152 
0153   cert_plot(threshold_pixel, threshold, filename, cert_strip, cert_track, cert_pixel, gLS_strip, gLS_track, gLS_pixel);
0154 }
0155 
0156 void ls_cert_type(string iDir,
0157                   float threshold,
0158                   string filename,
0159                   vector<string>& cert,
0160                   vector<pair<string, vector<float> > >& gLS,
0161                   vector<pair<string, vector<float> > >& bLS,
0162                   vector<pair<string, vector<float> > >& mLS) {
0163   void Cleaning(vector<int>&);
0164   string ListOut(vector<int>&);
0165 
0166   bool debug = false;
0167   string run = runnum_str(filename);
0168   if (debug)
0169     std::cout << filename.c_str() << std::endl;
0170 
0171   TDirectory* topDir;
0172   vector<float> ls;
0173 
0174   TFile* file = TFile::Open(filename.c_str());
0175   if (!file->IsOpen()) {
0176     std::cerr << "Failed to open " << filename << std::endl;
0177     return;
0178   }
0179 
0180   string dir = "DQMData/Run " + run + "/" + iDir;
0181   topDir = dynamic_cast<TDirectory*>(file->Get(dir.c_str()));
0182   topDir->cd();
0183   if (debug)
0184     std::cout << topDir->GetTitle() << std::endl;
0185 
0186   //
0187   // Reading the LS directory
0188   //
0189   TIter next(topDir->GetListOfKeys());
0190   TKey* key;
0191   while ((key = dynamic_cast<TKey*>(next()))) {
0192     string clName(key->GetClassName());
0193     if (clName == "TDirectoryFile") {
0194       TDirectory* curr_dir = dynamic_cast<TDirectory*>(key->ReadObj());
0195       string name = curr_dir->GetName();
0196       if (name == "Run summary")
0197         continue;
0198       name = name.substr(name.find('-') + 1);
0199       float temp1 = atof(name.c_str());
0200       ls.push_back(temp1);
0201     }
0202   }
0203   sort(ls.begin(), ls.end());
0204   int vecsize = ls.size();
0205 
0206   //
0207   // Definition of vectors for LS certification
0208   //
0209   Float_t* lsd = new Float_t[vecsize];
0210 
0211   Float_t** v = new Float_t*[4];
0212   for (int k = 0; k < 4; k++) {
0213     v[k] = new Float_t[vecsize];
0214   }
0215   //string certflag[4] = {"CertificationSummary","DAQSummary","DCSSummary","reportSummary"};
0216   //string certflagPrint[4] = {"Certification","DAQ","DCS","DQM"};
0217   string certflag[4] = {"DAQSummary", "DCSSummary", "reportSummary", "CertificationSummary"};
0218   string certflagPrint[4] = {"DAQ", "DCS", "DQM", "Certification"};
0219 
0220   int smax = 2;
0221   if (iDir == "SiStrip" || iDir == "Pixel")
0222     smax = 4;
0223 
0224   if (iDir == "Tracking") {
0225     certflag[0] = "CertificationSummary";
0226     certflagPrint[0] = "Certification";
0227     certflag[1] = "reportSummary";
0228     certflagPrint[1] = "DQM";
0229   }
0230 
0231   for (int icert_type = 0; icert_type < smax; icert_type++) {
0232     cert.push_back(certflagPrint[icert_type]);
0233   }
0234 
0235   if (debug)
0236     std::cout << gDirectory->GetName() << std::endl;
0237 
0238   for (int i = 0; i < vecsize; i++) {
0239     stringstream lsdir;
0240     lsdir << dir << "/By Lumi Section " << ls[i] << "-" << ls[i] << "/EventInfo";
0241     if (debug)
0242       std::cout << lsdir.str().c_str() << std::endl;
0243     float templs = ls[i];
0244     lsd[i] = templs;
0245     TDirectory* tempDir = dynamic_cast<TDirectory*>(file->Get(lsdir.str().c_str()));
0246     tempDir->cd();
0247     TIter nextTemp(tempDir->GetListOfKeys());
0248     TKey* keyTemp;
0249     while ((keyTemp = dynamic_cast<TKey*>(nextTemp()))) {
0250       float tempvalue = -1.;
0251       string classname(keyTemp->GetClassName());
0252       if (classname == "TObjString") {
0253         string sflag = keyTemp->GetName();
0254         string tempname = sflag.substr(sflag.find("f=") + 2);
0255         size_t pos1 = tempname.find('<');
0256         size_t pos2 = sflag.find_first_of('>');
0257         string detvalue = tempname.substr(0, pos1);
0258         string typecert = sflag.substr(1, pos2 - 1);
0259         if (debug)
0260           std::cout << typecert.c_str() << std::endl;
0261         tempvalue = atof(detvalue.c_str());
0262 
0263         for (int j = 0; j < smax; j++) {
0264           if (strstr(typecert.c_str(), certflag[j].c_str()) != nullptr)
0265             v[j][i] = tempvalue;
0266           if (debug)
0267             std::cout << "Entering value " << tempvalue << " " << v[j][i] << " for " << certflag[j].c_str()
0268                       << std::endl;
0269         }
0270       }
0271     }
0272   }
0273 
0274   int nLS_run = nlumis(filename);
0275 
0276   for (int iS = 0; iS < smax; iS++) {
0277     vector<int> goodLS;
0278     vector<int> badLS;
0279     vector<int> missingLS;
0280     vector<float> allLSthr;
0281 
0282     //loop over all available lumisections and fill good/bad lists
0283     for (int iLS = 0; iLS < vecsize; iLS++) {
0284       if (v[iS][iLS] >= threshold)
0285         goodLS.push_back(lsd[iLS]);
0286       else if (v[iS][iLS] > -1)  //protect from flagging non-tested LS as bad
0287         badLS.push_back(lsd[iLS]);
0288     }
0289 
0290     int last_idx = 0;
0291     for (int i_ls = 1; i_ls <= nLS_run; i_ls++) {
0292       for (int j = last_idx; j < vecsize; j++) {
0293         if (lsd[j] == i_ls) {
0294           last_idx = j + 1;
0295           if (v[iS][j] == 0)
0296             allLSthr.push_back(0.00001);
0297           else
0298             allLSthr.push_back(v[iS][j]);
0299           break;
0300         }
0301         if (lsd[j] > i_ls) {
0302           last_idx = j;
0303           missingLS.push_back(i_ls);
0304           allLSthr.push_back(-1);
0305           break;
0306         }
0307       }
0308     }
0309 
0310     Cleaning(goodLS);
0311     Cleaning(badLS);
0312     Cleaning(missingLS);
0313 
0314     string goodList = ListOut(goodLS);
0315     string badList = ListOut(badLS);
0316     string missingList = ListOut(missingLS);
0317 
0318     //save lumisections for this certification type
0319     gLS.push_back(make_pair(goodList, allLSthr));
0320     bLS.push_back(make_pair(badList, allLSthr));
0321     mLS.push_back(make_pair(missingList, allLSthr));
0322   }
0323 }
0324 
0325 void cert_plot(float threshold_pixel,
0326                float threshold,
0327                string filename,
0328                vector<string>& cert_strip,
0329                vector<string>& cert_track,
0330                vector<string>& cert_pixel,
0331                vector<pair<string, vector<float> > >& LS_strip,
0332                vector<pair<string, vector<float> > >& LS_track,
0333                vector<pair<string, vector<float> > >& LS_pixel) {
0334   int nLumiSections = nlumis(filename);
0335 
0336   char plottitles[200];
0337   sprintf(plottitles, "Lumisection Certification: Run %s;Luminosity Section;", runnum_str(filename).c_str());
0338   TH2D* cert_plot = new TH2D("cert_plot", plottitles, nLumiSections, 1, nLumiSections + 1, 5, 1, 6);
0339   cert_plot->SetStats(false);
0340   char label[100];
0341   for (int ityp = 0; ityp < 4; ityp++) {
0342     sprintf(label, "SiStrip %s", cert_strip[ityp].c_str());
0343     cert_plot->GetYaxis()->SetBinLabel(5 - ityp, label);
0344 
0345     for (unsigned int idx = 0; idx < LS_strip[ityp].second.size(); idx++)
0346       if (LS_strip[ityp].second[idx] > -1)
0347         cert_plot->SetBinContent(idx + 1, 5 - ityp, LS_strip[ityp].second[idx]);
0348   }
0349   for (int ityp = 0; ityp < 1; ityp++) {
0350     sprintf(label, "Tracking %s", cert_track[ityp].c_str());
0351     cert_plot->GetYaxis()->SetBinLabel(1 - ityp, label);
0352     for (unsigned int idx = 0; idx < LS_track[ityp].second.size(); idx++)
0353       if (LS_track[ityp].second[idx] > -1)
0354         cert_plot->SetBinContent(idx + 1, 1 - ityp, LS_track[ityp].second[idx]);
0355   }
0356 
0357   const Int_t colNum = 20;  // defining of a new palette
0358   Int_t palette[colNum];
0359   float rgb[colNum][3];
0360   int col_thr = colNum * threshold;
0361   for (Int_t i = 0; i < colNum; i++) {
0362     if (i >= col_thr) {
0363       // green
0364       rgb[i][0] = 0.00;
0365       rgb[i][1] = 0.80;
0366       rgb[i][2] = 0.00;
0367     } else {
0368       // red to yellow   //yellow red
0369       rgb[i][0] = 0.80 + (0.98 - 0.80) / (col_thr - 1) * i;  //0.98   0.80
0370       rgb[i][1] = 0.00 + (0.79 - 0.00) / (col_thr - 1) * i;  //0.79   0.00
0371       rgb[i][2] = 0.00;                                      //0.00
0372     }
0373 
0374     palette[i] = 9001 + i;
0375 
0376     TColor* color = gROOT->GetColor(9001 + i);
0377     if (!color)
0378       color = new TColor(9001 + i, 0, 0, 0, "");
0379     color->SetRGB(rgb[i][0], rgb[i][1], rgb[i][2]);
0380   }
0381   gStyle->SetPalette(colNum, palette);
0382   gROOT->SetStyle("Plain");
0383 
0384   TCanvas* cc = new TCanvas("name", "title", 1000, 600);
0385   cert_plot->Draw("colz");
0386   gPad->SetLeftMargin(0.17);
0387   string plotfilename = "Certification_run_" + runnum_str(filename) + ".png";
0388   cc->Print(plotfilename.c_str());
0389 
0390   //PIXEL plot
0391 
0392   TH2D* cert_plot_pixel = new TH2D("cert_plot_pixel", plottitles, nLumiSections, 1, nLumiSections + 1, 4, 1, 5);
0393   cert_plot_pixel->SetStats(false);
0394 
0395   for (int ityp = 0; ityp < 4; ityp++) {
0396     sprintf(label, "Pixel %s", cert_pixel[ityp].c_str());
0397     cert_plot_pixel->GetYaxis()->SetBinLabel(4 - ityp, label);
0398 
0399     for (unsigned int idx = 0; idx < LS_pixel[ityp].second.size(); idx++)
0400       if (LS_pixel[ityp].second[idx] > -1)
0401         cert_plot_pixel->SetBinContent(idx + 1, 4 - ityp, LS_pixel[ityp].second[idx]);
0402   }
0403 
0404   int col_thr_pixel = colNum * threshold_pixel;
0405 
0406   for (Int_t i = 0; i < colNum; i++) {
0407     if (i >= col_thr_pixel) {
0408       // green
0409       rgb[i][0] = 0.00;
0410       rgb[i][1] = 0.80;
0411       rgb[i][2] = 0.00;
0412     } else {
0413       // red to yellow   //yellow red
0414       rgb[i][0] = 0.80 + (0.98 - 0.80) / (col_thr - 1) * i;  //0.98   0.80
0415       rgb[i][1] = 0.00 + (0.79 - 0.00) / (col_thr - 1) * i;  //0.79   0.00
0416       rgb[i][2] = 0.00;                                      //0.00
0417     }
0418 
0419     palette[i] = 10001 + i;
0420 
0421     TColor* color = gROOT->GetColor(10001 + i);
0422     if (!color)
0423       color = new TColor(10001 + i, 0, 0, 0, "");
0424     color->SetRGB(rgb[i][0], rgb[i][1], rgb[i][2]);
0425   }
0426   gStyle->SetPalette(colNum, palette);
0427   gROOT->SetStyle("Plain");
0428 
0429   cert_plot_pixel->Draw("colz");
0430   //gPad->SetLeftMargin(0.17);
0431   string plotfilename_pixel = "Certification_run_" + runnum_str(filename) + "_pixel.png";
0432   cc->Print(plotfilename_pixel.c_str());
0433 
0434   delete cc;
0435 }
0436 
0437 int nlumis(string filename) {
0438   if (numlumis > -1)
0439     return numlumis;
0440 
0441   //TDirectory* topDir;
0442   vector<float> ls;
0443 
0444   TFile* file = TFile::Open(filename.c_str());
0445   if (!file->IsOpen()) {
0446     std::cerr << "Failed to open " << filename << std::endl;
0447     return -1;
0448   }
0449 
0450   string run = runnum_str(filename);
0451 
0452   //check if HIRun or pp run
0453   bool isHIRun = false;
0454   if (filename.find("HIRun") != string::npos)
0455     isHIRun = true;
0456 
0457   //valid up to the end of 2011 pp collisions
0458   if (!isHIRun) {
0459     string EventInfoDir = "DQMData/Run " + run + "/SiStrip/Run summary/EventInfo";
0460     TDirectory* rsEventInfoDir = dynamic_cast<TDirectory*>(file->Get(EventInfoDir.c_str()));
0461     rsEventInfoDir->cd();
0462     TIter eiKeys(rsEventInfoDir->GetListOfKeys());
0463     TKey* eiKey;
0464     while ((eiKey = dynamic_cast<TKey*>(eiKeys()))) {
0465       string classname(eiKey->GetClassName());
0466       if (classname == "TObjString") {
0467         string sflag = eiKey->GetName();
0468         string tempname = sflag.substr(sflag.find("i=") + 2);
0469         size_t pos1 = tempname.find('<');
0470         size_t pos2 = sflag.find_first_of('>');
0471         string detvalue = tempname.substr(0, pos1);
0472         string numlumisec = sflag.substr(1, pos2 - 1);
0473         if (numlumisec == (string) "iLumiSection") {
0474           numlumis = atoi(detvalue.c_str());
0475           break;
0476         }
0477       }
0478     }
0479   } else {
0480     //valid since 2011 HI running (iLumiSection variable not there anymore)
0481     string EventInfoDirHist = "DQMData/Run " + run + "/Info/Run summary/EventInfo/ProcessedLS";
0482     TH1F* allLS = (TH1F*)file->Get(EventInfoDirHist.c_str());
0483     numlumis = allLS->GetEntries() - 1;
0484     delete allLS;
0485   }
0486 
0487   return numlumis;
0488 }
0489 
0490 string runnum_str(string filename) { return filename.substr(filename.find("_R000") + 5, 6); }
0491 
0492 void Cleaning(vector<int>& LSlist) {
0493   if (LSlist.empty())
0494     return;
0495 
0496   //cleaning: keep only 1st and last lumisection in the range
0497   int refLS = LSlist[0];
0498   for (unsigned int at = 1; at < LSlist.size() - 1; at++) {
0499     //delete LSnums in between a single continuous range
0500     if (refLS + 1 == LSlist[at] && LSlist[at] + 1 == LSlist[at + 1]) {
0501       refLS = LSlist[at];
0502       LSlist[at] = -1;
0503     } else {
0504       refLS = LSlist[at];
0505     }
0506   }
0507 }
0508 
0509 string ListOut(vector<int>& LSlist) {
0510   string strout = "";
0511   bool rangeset = false;
0512 
0513   for (unsigned int at = 0; at < LSlist.size(); at++) {
0514     if (LSlist[at] != -1) {
0515       if (at > 0 && LSlist[at - 1] != -1)
0516         strout += ",";
0517       stringstream lsnum;
0518       lsnum << LSlist[at];
0519       strout += lsnum.str();
0520       rangeset = false;
0521     }
0522     if (LSlist[at] == -1 && !rangeset) {
0523       strout += "-";
0524       rangeset = true;
0525     }
0526   }
0527 
0528   return strout;
0529 }