ls_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
int numlumis = -1;

int     nlumis     ( string filename ); //get number of run lumisections
string  runnum_str ( string filename ); //read the run number, return in string
void    ls_cert( float threshold_pixel , float threshold , string filename ) ;

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

  if(argc==4) {
    char* cpixel = argv[1];
    char* cthr = argv[2];
    char* filename = argv[3];

    float threshold_pixel = 0;
    sscanf(cpixel,"%f",&threshold_pixel);
    float threshold = 0;
    sscanf(cthr,"%f",&threshold);

    std::cout << "ready to run ls_cert: pixel thr  " << threshold_pixel 
          << " threshold " << threshold
          << " filename " << filename << std::endl;

    ls_cert(threshold_pixel,threshold,filename);

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

}

void    ls_cert( float threshold_pixel , float threshold , string filename ) 
{
  void ls_cert_type ( string iDir , float threshold , string filename , vector <string>& , vector<pair<string,vector<float> > >& , vector<pair<string,vector<float> > >& , vector<pair<string,vector<float> > >& );
  void cert_plot    ( float threshold_pixel , float threshold , string filename , 
              vector <string>& , vector <string>& , vector <string>& , vector<pair<string,vector<float> > >& , vector<pair<string,vector<float> > >& , vector<pair<string,vector<float> > >& );

  //presets
  numlumis = -1;

  //certifications
  vector <string> cert_strip;
  vector <string> cert_track;
  vector <string> cert_pixel;
  
  //good lumisections
  vector < pair < string,vector <float> > > gLS_strip;
  vector < pair < string,vector <float> > > gLS_track;
  vector < pair < string,vector <float> > > gLS_pixel;

  //bad lumisections
  vector < pair < string,vector <float> > > bLS_strip;
  vector < pair < string,vector <float> > > bLS_track;
  vector < pair < string,vector <float> > > bLS_pixel;

  //missing lumisections
  vector < pair < string,vector <float> > > mLS_strip;
  vector < pair < string,vector <float> > > mLS_track;
  vector < pair < string,vector <float> > > mLS_pixel;

  ls_cert_type( "SiStrip"  , threshold       , filename , cert_strip , gLS_strip , bLS_strip , mLS_strip );
  ls_cert_type( "Tracking" , threshold       , filename , cert_track , gLS_track , bLS_track , mLS_track );
  ls_cert_type( "Pixel"    , threshold_pixel , filename , cert_pixel , gLS_pixel , bLS_pixel , mLS_pixel );

  std::ofstream outfile;
  string namefile = "Certification_run_" + runnum_str( filename ) + ".txt";
  outfile.open(namefile.c_str());
  outfile << "Lumisection Certification (GOOD: >= " << threshold_pixel << " [Pixel]; >= " << threshold << " [SiStrip,Tracking] " 
      << ", otherwise BAD):" << endl << endl;
  outfile << "GOOD Lumisections:" << endl;
  char line[200];
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " Pixel    %*sSummary: %s" , 13 , cert_pixel[ityp].c_str() , gLS_pixel[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " SiStrip  %*sSummary: %s" , 13 , cert_strip[ityp].c_str() , gLS_strip[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 1; ityp++)
    {
      sprintf( line, " Tracking %*sSummary: %s" , 13 , cert_track[ityp].c_str() , gLS_track[ityp].first.c_str() );
      outfile << line << endl;
    }

  outfile << "\nBAD Lumisections:" << endl;
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " Pixel    %*sSummary: %s" , 13 , cert_pixel[ityp].c_str() , bLS_pixel[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " SiStrip  %*sSummary: %s" , 13 , cert_strip[ityp].c_str() , bLS_strip[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 1; ityp++)
    {
      sprintf( line, " Tracking %*sSummary: %s" , 13 , cert_track[ityp].c_str() , bLS_track[ityp].first.c_str() );
      outfile << line << endl;
    }
  
  outfile << "\nMISSING Lumisections:" << endl;
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " Pixel    %*sSummary: %s" , 13 , cert_pixel[ityp].c_str() , mLS_pixel[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( line, " SiStrip  %*sSummary: %s" , 13 , cert_strip[ityp].c_str() , mLS_strip[ityp].first.c_str() );
      outfile << line << endl;
    }
  for (int ityp = 0; ityp < 1; ityp++)
    {
      sprintf( line, " Tracking %*sSummary: %s" , 13 , cert_track[ityp].c_str() , mLS_track[ityp].first.c_str() );
      outfile << line << endl;
    }
  
  outfile.close();
  std::cout << "Lumisection Certification summary saved in " << namefile << std::endl;

  cert_plot ( threshold_pixel , threshold , filename , cert_strip , cert_track , cert_pixel , gLS_strip , gLS_track , gLS_pixel );
}

void ls_cert_type(string iDir, float threshold, string filename, vector <string>& cert, vector<pair<string,vector<float> > >& gLS, vector<pair<string,vector<float> > >& bLS, vector<pair<string,vector<float> > >& mLS ) 
{  
  void Cleaning(vector<int> &);
  string ListOut(vector<int> &);

  bool debug = false;
  string run = runnum_str( filename );
  if (debug) std::cout << filename.c_str() << std::endl;

  TDirectory* topDir; 
  vector<float> ls;
  
  TFile* file = TFile::Open(filename.c_str());
  if (!file->IsOpen()) {
    std::cerr << "Failed to open " << filename << std::endl; 
    return;
  }

  string dir = "DQMData/Run " + run + "/" + iDir;
  topDir = dynamic_cast<TDirectory*>( file->Get(dir.c_str()));
  topDir->cd();
  if (debug) std::cout << topDir->GetTitle() << std::endl;

  //
  // Reading the LS directory    
  //
  TIter next(topDir->GetListOfKeys());
  TKey *key;
  while  ( (key = dynamic_cast<TKey*>(next())) ) {
    string clName(key->GetClassName());
    if (clName == "TDirectoryFile") {
      TDirectory *curr_dir = dynamic_cast<TDirectory*>(key->ReadObj());
      string name = curr_dir->GetName();
      if (name == "Run summary") continue;
      name = name.substr(name.find("-")+1);
      float temp1 = atof(name.c_str()); 
      ls.push_back(temp1);
    }
  }
  sort(ls.begin(),ls.end());   
  int vecsize = ls.size();

  //
  // Definition of vectors for LS certification
  //
  Float_t * lsd = new Float_t[vecsize];
  
  Float_t ** v = new Float_t*[4];
  for(int k=0;k<4;k++) {
    v[k] = new Float_t[vecsize];
  }
  //string certflag[4] = {"CertificationSummary","DAQSummary","DCSSummary","reportSummary"};
  //string certflagPrint[4] = {"Certification","DAQ","DCS","DQM"};
  string certflag[4] = {"DAQSummary","DCSSummary","reportSummary","CertificationSummary"};
  string certflagPrint[4] = {"DAQ","DCS","DQM","Certification"};

  int smax = 2;
  if ( iDir == "SiStrip" || iDir == "Pixel" ) smax = 4;
  
  if ( iDir == "Tracking" )
    {
      certflag[0] = "CertificationSummary";
      certflagPrint[0] = "Certification";
      certflag[1] = "reportSummary";
      certflagPrint[1] = "DQM";
    }

  for (int icert_type = 0; icert_type < smax; icert_type++)
    {
      cert.push_back( certflagPrint[icert_type] );
    }

  if (debug) std::cout << gDirectory->GetName() << std::endl;
  
  for (int i=0; i < vecsize; i++){
    stringstream lsdir;
    lsdir << dir << "/By Lumi Section " << ls[i] <<"-"<<ls[i]<<"/EventInfo";
    if (debug) std::cout << lsdir.str().c_str() << std::endl;
    float templs = ls[i];
    lsd[i] = templs;
    TDirectory *tempDir = dynamic_cast<TDirectory*>( file->Get(lsdir.str().c_str()));
    tempDir->cd();
    int j = 0;      
    TIter nextTemp(tempDir->GetListOfKeys());
    TKey *keyTemp;
    while  ( (keyTemp = dynamic_cast<TKey*>(nextTemp())) ) {
      float tempvalue = -1.;
      string classname(keyTemp->GetClassName());
      if (classname=="TObjString" ){
    string sflag = keyTemp->GetName();
    string tempname = sflag.substr(sflag.find("f=")+2);
    size_t pos1 = tempname.find("<");
    size_t pos2 = sflag.find_first_of(">");
    string detvalue = tempname.substr(0,pos1);
    string typecert = sflag.substr(1,pos2-1);
    if (debug) std::cout << typecert.c_str() << std::endl;
    tempvalue = atof(detvalue.c_str());
    
    for (j=0; j<smax; j++){
      if ( strstr(typecert.c_str(),certflag[j].c_str())!=nullptr)
        v[j][i] = tempvalue;
      if (debug) std::cout << "Entering value " << tempvalue << " " << v[j][i] << " for " << certflag[j].c_str() << std::endl;
    }
    j = j + 1;
      }
    }
  }

  int nLS_run = nlumis ( filename );

  for (int iS = 0; iS < smax; iS++)
    {
      
      vector<int> goodLS;
      vector<int> badLS;
      vector<int> missingLS;
      vector<float> allLSthr;

      //loop over all available lumisections and fill good/bad lists
      for (int iLS = 0; iLS < vecsize; iLS++)
    {
      if ( v[iS][iLS] >= threshold ) 
        goodLS.push_back(lsd[iLS]);
      else                  
        if ( v[iS][iLS] > -1 ) //protect from flagging non-tested LS as bad
          badLS.push_back(lsd[iLS]);
    }

      int last_idx = 0;
      for (int i_ls = 1; i_ls <= nLS_run; i_ls++)
    {
      for (int j = last_idx; j < vecsize; j++)
        {
          if ( lsd[j] == i_ls ) 
        {
          last_idx = j+1;
          if ( v[iS][j] == 0 ) allLSthr.push_back(0.00001);
          else
            allLSthr.push_back(v[iS][j]);
          break;
        }
          if ( lsd[j] > i_ls )
        {
          last_idx = j;
          missingLS.push_back(i_ls);
          allLSthr.push_back(-1);
          break;
        }
        }
    }

      Cleaning( goodLS    );
      Cleaning( badLS     );
      Cleaning( missingLS );

      string goodList    = ListOut( goodLS    );
      string badList     = ListOut( badLS     );
      string missingList = ListOut( missingLS );

      //save lumisections for this certification type
      gLS.push_back ( make_pair ( goodList    , allLSthr ) );
      bLS.push_back ( make_pair ( badList     , allLSthr ) );
      mLS.push_back ( make_pair ( missingList , allLSthr ) );      
    }
}

void cert_plot( float threshold_pixel , float threshold , string filename , vector <string>& cert_strip , vector <string>& cert_track , 
        vector <string>& cert_pixel , vector<pair<string,vector<float> > >& LS_strip , vector<pair<string,vector<float> > >& LS_track , 
        vector<pair<string,vector<float> > >& LS_pixel )
{
  int nLumiSections = nlumis( filename );

  char plottitles[200];
  sprintf( plottitles , "Lumisection Certification: Run %s;Luminosity Section;" , runnum_str( filename ).c_str() );
  TH2D *cert_plot = new TH2D( "cert_plot" , plottitles , nLumiSections , 1 , nLumiSections + 1 , 5 , 1 , 6 );
  cert_plot->SetStats(false);
  char label[100];
  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( label , "SiStrip %s" , cert_strip[ityp].c_str() );
      cert_plot->GetYaxis()->SetBinLabel( 5 - ityp , label );

      for (unsigned int idx = 0; idx <  LS_strip[ityp].second.size() ; idx++)
        if ( LS_strip[ityp].second[idx] > -1 ) cert_plot->SetBinContent(idx+1, 5 - ityp, LS_strip[ityp].second[idx]);
    }
  for (int ityp = 0; ityp < 1; ityp++)
    {
      sprintf( label , "Tracking %s" , cert_track[ityp].c_str() );
      cert_plot->GetYaxis()->SetBinLabel( 1 - ityp , label );
      for (unsigned int idx = 0; idx <  LS_track[ityp].second.size() ; idx++)
        if ( LS_track[ityp].second[idx] > -1 ) cert_plot->SetBinContent(idx+1, 1 - ityp, LS_track[ityp].second[idx]);
    }

  const Int_t colNum = 20; // defining of a new palette
  Int_t palette[colNum];
  float rgb[colNum][3];
  int col_thr = colNum * threshold;
  for (Int_t i=0; i<colNum; i++)
    {
      if ( i >= col_thr )
        {
          // green
          rgb[i][0] = 0.00;
          rgb[i][1] = 0.80;
          rgb[i][2] = 0.00;
        }
      else
        {
          // red to yellow   //yellow red
          rgb[i][0] = 0.80 + ( 0.98 - 0.80 ) / ( col_thr - 1 ) * i ;  //0.98   0.80
          rgb[i][1] = 0.00 + ( 0.79 - 0.00 ) / ( col_thr - 1 ) * i ;  //0.79   0.00
          rgb[i][2] = 0.00;  //0.00
        }

      palette[i] = 9001+i;

      TColor *color = gROOT->GetColor(9001+i);
      if (!color) color = new TColor(9001 + i, 0, 0, 0, "");
      color->SetRGB(rgb[i][0], rgb[i][1], rgb[i][2]);
    }
  gStyle->SetPalette(colNum,palette);
  gROOT->SetStyle("Plain");

  TCanvas *cc = new TCanvas( "name" , "title" , 1000 , 600 );
  cert_plot->Draw("colz");
  gPad->SetLeftMargin(0.17);
  string plotfilename = "Certification_run_" +  runnum_str( filename ) + ".png";
  cc->Print( plotfilename.c_str() );


  //PIXEL plot

  TH2D *cert_plot_pixel = new TH2D( "cert_plot_pixel" , plottitles , nLumiSections , 1 , nLumiSections + 1 , 4 , 1 , 5 );
  cert_plot_pixel->SetStats(false);

  for (int ityp = 0; ityp < 4; ityp++)
    {
      sprintf( label , "Pixel %s" , cert_pixel[ityp].c_str() );
      cert_plot_pixel->GetYaxis()->SetBinLabel( 4 - ityp , label );

      for (unsigned int idx = 0; idx <  LS_pixel[ityp].second.size() ; idx++)
        if ( LS_pixel[ityp].second[idx] > -1 ) cert_plot_pixel->SetBinContent(idx+1, 4 - ityp, LS_pixel[ityp].second[idx]);
    }

  int col_thr_pixel = colNum * threshold_pixel;

  for (Int_t i=0; i<colNum; i++)
    {
      if ( i >= col_thr_pixel )
        {
          // green
          rgb[i][0] = 0.00;
          rgb[i][1] = 0.80;
          rgb[i][2] = 0.00;
        }
      else
        {
          // red to yellow   //yellow red
          rgb[i][0] = 0.80 + ( 0.98 - 0.80 ) / ( col_thr - 1 ) * i ;  //0.98   0.80
          rgb[i][1] = 0.00 + ( 0.79 - 0.00 ) / ( col_thr - 1 ) * i ;  //0.79   0.00
          rgb[i][2] = 0.00;  //0.00
        }

      palette[i] = 10001+i;

      TColor *color = gROOT->GetColor(10001+i);
      if (!color) color = new TColor(10001 + i, 0, 0, 0, "");
      color->SetRGB(rgb[i][0], rgb[i][1], rgb[i][2]);
    }
  gStyle->SetPalette(colNum,palette);
  gROOT->SetStyle("Plain");

  cert_plot_pixel->Draw("colz");
  //gPad->SetLeftMargin(0.17);
  string plotfilename_pixel = "Certification_run_" +  runnum_str( filename ) + "_pixel.png";
  cc->Print( plotfilename_pixel.c_str() );

  delete cc;
}


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

  //TDirectory* topDir;
  vector<float> ls;

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

  //check if HIRun or pp run
  bool isHIRun = false;
  if ( filename.find ( "HIRun" ) != string::npos )
    isHIRun = true;

  //valid up to the end of 2011 pp collisions
  if ( !isHIRun )
    {
      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;
        }
        }
    }
    }
  else
    {
      //valid since 2011 HI running (iLumiSection variable not there anymore)  
      string EventInfoDirHist = "DQMData/Run " + run + "/Info/Run summary/EventInfo/ProcessedLS";
      TH1F* allLS = (TH1F*) file->Get( EventInfoDirHist.c_str() );
      numlumis = allLS->GetEntries() - 1;
      delete allLS;
    }

  return numlumis;
}

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

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)
{
  
  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;
}