SiPixelInformationExtractor.cc

Go to the documentation of this file.

#include "DQM/SiPixelMonitorClient/interface/SiPixelInformationExtractor.h"
#include "DQM/SiPixelMonitorClient/interface/SiPixelUtility.h"
#include "DQM/SiPixelMonitorClient/interface/SiPixelEDAClient.h"
#include "DQM/SiPixelMonitorClient/interface/ANSIColors.h"
#include "DQM/SiPixelCommon/interface/SiPixelFolderOrganizer.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "DataFormats/SiPixelDetId/interface/PixelBarrelName.h"
#include "DataFormats/SiPixelDetId/interface/PixelBarrelNameUpgrade.h"
#include "DataFormats/SiPixelDetId/interface/PixelEndcapName.h"
#include "DataFormats/SiPixelDetId/interface/PixelEndcapNameUpgrade.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"

#include "CondFormats/SiPixelObjects/interface/DetectorIndex.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelFrameConverter.h"

#include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"

#include "TClass.h"
#include "TText.h"
#include "TROOT.h"
#include "TPad.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TString.h"
#include "TImage.h"
#include "TPaveText.h"
#include "TImageDump.h"
#include "TRandom.h"
#include "TStopwatch.h"
#include "TAxis.h"
#include "TPaveLabel.h"
#include "Rtypes.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"

#include <iostream>
#include <cmath>
#include <map>

#include <cstdlib> // for free() - Root can allocate with malloc() - sigh...
 
using namespace std;
using namespace edm;

//------------------------------------------------------------------------------
SiPixelInformationExtractor::SiPixelInformationExtractor(bool offlineXMLfile) : offlineXMLfile_(offlineXMLfile) {
  edm::LogInfo("SiPixelInformationExtractor") << 
    " Creating SiPixelInformationExtractor " << "\n" ;
  
  readReference_ = false;
}

//------------------------------------------------------------------------------
SiPixelInformationExtractor::~SiPixelInformationExtractor() {
  edm::LogInfo("SiPixelInformationExtractor") << 
    " Deleting SiPixelInformationExtractor " << "\n" ;
  
}

//------------------------------------------------------------------------------
void SiPixelInformationExtractor::readConfiguration() { }


//============================================================================================================
// --  Return type of ME
//
std::string  SiPixelInformationExtractor::getMEType(MonitorElement * theMe)
{
  string qtype = theMe->getRootObject()->IsA()->GetName() ;
  if(         qtype.find("TH1") != string::npos )
  {
    return "TH1" ;
  } else if ( qtype.find("TH2") != string::npos  ) {
    return "TH2" ;
  } else if ( qtype.find("TH3") != string::npos ) {
    return "TH3" ;
  }
  return "TH1" ;
}

//------------------------------------------------------------------------------
void SiPixelInformationExtractor::getItemList(const multimap<string, string>& req_map, 
                                              string item_name,
                          vector<string>& items) {
  items.clear();
  for (multimap<string, string>::const_iterator it = req_map.begin();
       it != req_map.end(); it++) {
    if (it->first == item_name) {
      items.push_back(it->second);
    }
  }
}

//------------------------------------------------------------------------------
bool SiPixelInformationExtractor::hasItem(multimap<string,string>& req_map,
                      string item_name){
  multimap<string,string>::iterator pos = req_map.find(item_name);
  if (pos != req_map.end()) return true;
  return false;  
}

//------------------------------------------------------------------------------
std::string SiPixelInformationExtractor::getItemValue(const std::multimap<std::string,std::string>& req_map,
                         std::string item_name){
  std::multimap<std::string,std::string>::const_iterator pos = req_map.find(item_name);
  std::string value = " ";
  if (pos != req_map.end()) {
    value = pos->second;
  }
  return value;
}
std::string SiPixelInformationExtractor::getItemValue(std::multimap<std::string,std::string>& req_map,
                         std::string item_name){
  std::multimap<std::string,std::string>::iterator pos = req_map.find(item_name);
  std::string value = " ";
  if (pos != req_map.end()) {
    value = pos->second;
  }
  return value;
}

//
// -- Get color  name from status
//
void SiPixelInformationExtractor::selectColor(string& col, int status){
  if (status == dqm::qstatus::STATUS_OK)    col = "#00ff00";
  else if (status == dqm::qstatus::WARNING) col = "#ffff00";
  else if (status == dqm::qstatus::ERROR)   col = "#ff0000";
  else if (status == dqm::qstatus::OTHER)   col = "#ffa500";
  else  col = "#0000ff";
}
//
// -- Get Image name from ME
//
void SiPixelInformationExtractor::selectColor(string& col, vector<QReport*>& reports){
  int istat = 999;
  int status = 0;
  for (vector<QReport*>::const_iterator it = reports.begin(); it != reports.end();
       it++) {
    status = (*it)->getStatus();
    if (status > istat) istat = status;
  }
  selectColor(col, status);
}
//
// -- Get Image name from status
//
void SiPixelInformationExtractor::selectImage(string& name, int status){
  if (status == dqm::qstatus::STATUS_OK) name="images/LI_green.gif";
  else if (status == dqm::qstatus::WARNING) name="images/LI_yellow.gif";
  else if (status == dqm::qstatus::ERROR) name="images/LI_red.gif";
  else if (status == dqm::qstatus::OTHER) name="images/LI_orange.gif";
  else  name="images/LI_blue.gif";
}
//
// -- Get Image name from ME
//
void SiPixelInformationExtractor::selectImage(string& name, vector<QReport*>& reports){
  int istat = 999;
  int status = 0;
  for (vector<QReport*>::const_iterator it = reports.begin(); it != reports.end();
       it++) {
    status = (*it)->getStatus();
    if (status > istat) istat = status;
  }
  selectImage(name, status);
}

//------------------------------------------------------------------------------
void SiPixelInformationExtractor::computeStatus(MonitorElement      * theME,
                                                double              & colorValue,
                        pair<double,double> & norm) 
{
  double normalizationX = 1 ;
  double normalizationY = 1 ;
  double meanX          = 0 ;
  double meanY          = 0 ;
  
  colorValue = 0 ;

  pair<double,double> normX ;
  pair<double,double> normY ;

  string theMEType = getMEType(theME) ;

  if( theMEType.find("TH1") != string::npos)
  {
   meanX = (double)theME->getMean();
   getNormalization(theME, normX, "TH1") ;
   normalizationX = fabs( normX.second - normX.first) ;
   if( normalizationX == 0 ) {normalizationX=1.E-20;}
   colorValue  = meanX / normalizationX ;
   norm.first  = normX.first ;
   norm.second = normX.second ;
  }
  
  if( theMEType.find("TH2") != string::npos)
  {
   meanX = (double)theME->getMean(1);
   meanY = (double)theME->getMean(2);
   getNormalization2D(theME, normX, normY, "TH2") ;
   normalizationX = fabs( normX.second - normX.first) ;
   normalizationY = fabs( normY.second - normY.first) ;
   if( normalizationX == 0 ) {normalizationX=1.E-20;}
   if( normalizationY == 0 ) {normalizationY=1.E-20;}
   double cVX = meanX / normalizationX ;
   double cVY = meanY / normalizationY ;
   colorValue = sqrt(cVX*cVX + cVY*cVY) ;
   if( normalizationX >= normalizationY )
   { 
    norm.first  = normX.first;
    norm.second = normX.second ;
   } else { 
    norm.first  = normY.first;
    norm.second = normY.second ;
   }
  } 
 
  return ;
}

//------------------------------------------------------------------------------
void SiPixelInformationExtractor::getNormalization(MonitorElement     * theME, 
                                                   pair<double,double>& norm,
                           std::string          theMEType) 
{
  double normLow  = 0 ;
  double normHigh = 0 ;

  if( theMEType.find("TH1") != string::npos)
  {
   normHigh    = (double)theME->getNbinsX() ;
   norm.first  = normLow  ;
   norm.second = normHigh ;
  }
}

//------------------------------------------------------------------------------
void SiPixelInformationExtractor::getNormalization2D(MonitorElement     * theME, 
                                                     pair<double,double>& normX,
                                                     pair<double,double>& normY,
                             std::string          theMEType) 
{
  double normLow  = 0 ;
  double normHigh = 0 ;

  if( theMEType.find("TH2") != string::npos )
  {
   normHigh    = (double)theME->getNbinsX() ;
   normX.first  = normLow  ;
   normX.second = normHigh ;
   normHigh    = (double)theME->getNbinsY() ;
   normY.first  = normLow  ;
   normY.second = normHigh ;
  }
}

//------------------------------------------------------------------------------
int SiPixelInformationExtractor::getDetId(MonitorElement * mE) 
{
 const string& mEName = mE->getName();

 int detId = 0;
 
 if( mEName.find("_3") != string::npos )
 {
  string detIdString = mEName.substr((mEName.find_last_of("_"))+1,9);
  std::istringstream isst;
  isst.str(detIdString);
  isst>>detId;
 }
  return detId ;
  
}


void SiPixelInformationExtractor::bookNoisyPixels(DQMStore::IBooker & iBooker, float noiseRate_,bool Tier0Flag) {
//std::cout<<"BOOK NOISY PIXEL MEs!"<<std::endl;
  iBooker.cd();
  if(noiseRate_>=0.){
    iBooker.setCurrentFolder("Pixel/Barrel");
    EventRateBarrelPixels = iBooker.book1D("barrelEventRate","Digi event rate for all Barrel pixels",1000,0.,0.01);
    EventRateBarrelPixels->setAxisTitle("Event Rate",1);
    EventRateBarrelPixels->setAxisTitle("Number of Pixels",2);
    iBooker.cd();  
    iBooker.setCurrentFolder("Pixel/Endcap");
    EventRateEndcapPixels = iBooker.book1D("endcapEventRate","Digi event rate for all Endcap pixels",1000,0.,0.01);
    EventRateEndcapPixels->setAxisTitle("Event Rate",1);
    EventRateEndcapPixels->setAxisTitle("Number of Pixels",2);
  }
}

          

void SiPixelInformationExtractor::findNoisyPixels(DQMStore::IBooker & iBooker, DQMStore::IGetter & iGetter, bool init, float noiseRate_, int noiseRateDenominator_, edm::ESHandle<SiPixelFedCablingMap> theCablingMap)
{

  
  if(init){
    endOfModules_=false;
    nevents_=noiseRateDenominator_;
    if(nevents_ == -1){
      iBooker.cd();
      iGetter.cd();
      iBooker.setCurrentFolder("Pixel/EventInfo");
      iGetter.setCurrentFolder("Pixel/EventInfo");
      nevents_ = (iGetter.get("Pixel/EventInfo/processedEvents"))->getIntValue();
    }
    iBooker.cd();
    iGetter.cd();
    myfile_.open ("NoisyPixelList.txt", ios::app);
    myfile_ << "Noise summary, ran over " << nevents_ << " events, threshold was set to " << noiseRate_ <<  std::endl;
  }
  string currDir = iBooker.pwd();
  string dname = currDir.substr(currDir.find_last_of("/")+1);


  if(dname.find("Module_")!=string::npos){
    vector<string> meVec = iGetter.getMEs();
    for (vector<string>::const_iterator it = meVec.begin(); it != meVec.end(); it++) {
      string full_path = currDir + "/" + (*it);
      if(full_path.find("hitmap_siPixelDigis")!=string::npos){
        MonitorElement * me = iGetter.get(full_path);
        if (!me) continue;
    int detid=getDetId(me); int pixcol=-1; int pixrow=-1; 
    std::vector<std::pair<std::pair<int, int>, float> > noisyPixelsInModule;
    TH2F * hothisto = me->getTH2F();
    if(hothisto){
      for(int i=1; i!=hothisto->GetNbinsX()+1; i++){
        for(int j=1; j!=hothisto->GetNbinsY()+1; j++){
          float value = (hothisto->GetBinContent(i,j))/float(nevents_);
          if(me->getPathname().find("Barrel")!=string::npos){
            EventRateBarrelPixels = iGetter.get("Pixel/Barrel/barrelEventRate");
            if(EventRateBarrelPixels) EventRateBarrelPixels->Fill(value);
          }else if(me->getPathname().find("Endcap")!=string::npos){
            EventRateEndcapPixels = iGetter.get("Pixel/Endcap/endcapEventRate");
            if(EventRateEndcapPixels) EventRateEndcapPixels->Fill(value);
          }
          if(value > noiseRate_){
            pixcol = i-1;
            pixrow = j-1;
 
            std::pair<int, int> address(pixcol, pixrow);
            std::pair<std::pair<int, int>, float>  PixelStats(address, value);
            noisyPixelsInModule.push_back(PixelStats);
          }
            }
      }
    }
    noisyDetIds_[detid] = noisyPixelsInModule;
      }
    }
  }
  vector<string> subDirVec = iGetter.getSubdirs();  
  for (vector<string>::const_iterator ic = subDirVec.begin();
       ic != subDirVec.end(); ic++) {
    if((*ic).find("AdditionalPixelErrors")!=string::npos) continue;
    iGetter.cd(*ic);
    iBooker.cd(*ic);
    init=false;
    findNoisyPixels(iBooker,iGetter,init,noiseRate_,noiseRateDenominator_,theCablingMap);
    iBooker.goUp();
    iGetter.setCurrentFolder(iBooker.pwd());
  }

  if(iBooker.pwd().find("EventInfo")!=string::npos) endOfModules_ = true;
  
  if(!endOfModules_) return;
  if(currDir == "Pixel/EventInfo/reportSummaryContents"){
    std::vector<std::pair<sipixelobjects::DetectorIndex,double> > pixelvec;
    std::map<uint32_t,int> myfedmap;
    std::map<uint32_t,std::string> mynamemap;
    int realfedID = -1;
    int counter = 0;
    int n_noisyrocs_all = 0;
    int n_noisyrocs_barrel = 0;
    int n_noisyrocs_endcap = 0;
    int n_verynoisyrocs_all = 0;
    int n_verynoisyrocs_barrel = 0;
    int n_verynoisyrocs_endcap = 0;

    for(int fid = 0; fid < 40; fid++){
    for(std::map<uint32_t, std::vector< std::pair<std::pair<int, int>, float> > >::const_iterator it = noisyDetIds_.begin(); 
        it != noisyDetIds_.end(); it++){
      uint32_t detid = (*it).first;
      std::vector< std::pair<std::pair<int, int>, float> > noisyPixels = (*it).second;
      // now convert into online conventions:
      for(int fedid=0; fedid<=40; ++fedid){
    SiPixelFrameConverter converter(theCablingMap.product(),fedid);
    uint32_t newDetId = detid;
    if(converter.hasDetUnit(newDetId)){
      realfedID=fedid;
      break;   
    }
      }
      if(fid == realfedID){
      if(realfedID==-1) continue; 
      DetId detId(detid);
      uint32_t detSubId = detId.subdetId();
      std::string outputname;
      bool HalfModule = false;
      if (detSubId == 2){   //FPIX
    PixelEndcapName nameworker(detid);
    outputname = nameworker.name();
      } else if(detSubId == 1){   //BPIX
    PixelBarrelName nameworker(detid);
    outputname = nameworker.name();
    HalfModule = nameworker.isHalfModule();

      } else{
    continue;
      } 
      std::map<int,int> myrocmap;
      myfedmap[detid]=realfedID;
      mynamemap[detid]=outputname;
      
      for(std::vector< std::pair< std::pair<int,int>, float> >::const_iterator pxl = noisyPixels.begin(); 
          pxl != noisyPixels.end(); pxl++){
        std::pair<int,int> offlineaddress = (*pxl).first;
    float Noise_frac = (*pxl).second;
    int offlineColumn = offlineaddress.first;
        int offlineRow = offlineaddress.second;
        counter++;

        sipixelobjects::ElectronicIndex cabling; 
        SiPixelFrameConverter formatter(theCablingMap.product(),realfedID);
        sipixelobjects::DetectorIndex detector = {detid, offlineRow, offlineColumn};      
    formatter.toCabling(cabling,detector);
        // cabling should now contain cabling.roc and cabling.dcol  and cabling.pxid
        // however, the coordinates now need to be converted from dcl,pxid to the row,col coordinates used in the calibration info 
        sipixelobjects::LocalPixel::DcolPxid loc;
        loc.dcol = cabling.dcol;
        loc.pxid = cabling.pxid;
    
            sipixelobjects::LocalPixel locpixel(loc);
        assert(realfedID >= 0);
        assert(cabling.link >= 0);
        assert(cabling.roc >= 0);
    sipixelobjects::CablingPathToDetUnit path = {static_cast<unsigned int>(realfedID), 
                                                     static_cast<unsigned int>(cabling.link),
                                                     static_cast<unsigned int>(cabling.roc)};  
    const sipixelobjects::PixelROC *theRoc = theCablingMap->findItem(path);
    // END of FIX
    
        int onlineColumn = locpixel.rocCol();
        int onlineRow= locpixel.rocRow();
    myrocmap[(theRoc->idInDetUnit())]++;

    // ROC numbers in the barrel go from 8 to 15 instead of 0 to 7 in half modules.  This is a 
    // fix to get the roc number, and add 8 to it if:
    // it's a Barrel module AND on the minus side AND a Half module

    int rocnumber = -1;

    if((detSubId == 1) && (outputname.find("mO")!=string::npos || outputname.find("mI")!=string::npos) && (HalfModule)){
      rocnumber = theRoc->idInDetUnit() + 8;
    }
    else{
      rocnumber = theRoc->idInDetUnit();
    }

        myfile_ <<"NAME: "<<outputname<<" , DETID: "<<detid<<" , OFFLINE: col,row: "<<offlineColumn<<","<<offlineRow<<"  \t , ONLINE: roc,col,row: "<<rocnumber<<","<<onlineColumn<<","<<onlineRow<< "  \t , fed,dcol,pixid,link: "<<realfedID<<","<<loc.dcol<<","<<loc.pxid<<","<<cabling.link << ", Noise fraction: " << Noise_frac << std::endl;
      }
      for(std::map<int, int>::const_iterator nrc = myrocmap.begin(); nrc != myrocmap.end(); nrc++){
    if((*nrc).second > 0){
      n_noisyrocs_all++;
      if(detSubId == 2){
        n_noisyrocs_endcap++;
      } else if(detSubId == 1){
        n_noisyrocs_barrel++;}
    }
    if((*nrc).second > 40){
      n_verynoisyrocs_all++;
      if(detSubId == 2){
        n_verynoisyrocs_endcap++;
      } else if(detSubId == 1){
        n_verynoisyrocs_barrel++;}
    }
      }
      }

    }
    }
    myfile_ << "There are " << n_noisyrocs_all << " noisy ROCs (ROCs with at least 1 noisy pixel) in the entire detector. " << n_noisyrocs_endcap << " are in the FPIX and " << n_noisyrocs_barrel << " are in the BPIX. " << endl;
    myfile_ << "There are " << n_verynoisyrocs_all << " highly noisy ROCs (ROCs with at least 10% of all pixels passing the noise threshold) in the entire detector. " << n_verynoisyrocs_endcap << " are in the FPIX and " << n_verynoisyrocs_barrel << " are in the BPIX. " << endl;

  }
  myfile_.close();
  return;
}