EmDQMPostProcessor.cc

Go to the documentation of this file.

#include "HLTriggerOffline/Egamma/interface/EmDQMPostProcessor.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Framework/interface/MakerMacros.h"


#include <iostream>
#include <string.h>
#include <iomanip>
#include <fstream>
#include <math.h>


EmDQMPostProcessor::EmDQMPostProcessor(const edm::ParameterSet& pset)
{
  subDir_ = pset.getUntrackedParameter<std::string>("subDir");
}

void EmDQMPostProcessor::endRun(edm::Run const& run, edm::EventSetup const& es)
{
  // setup DQM stor               //
  
  DQMStore * dqm = 0;
  dqm = edm::Service<DQMStore>().operator->();

  if ( ! dqm ) {
    edm::LogInfo("EmDQMPostProcessor") << "Cannot create DQMStore instance\n";
    return;
  }


  //go to the directory to be processed
  if(dqm->dirExists(subDir_)) dqm->cd(subDir_);
  else {
    edm::LogWarning("EmDQMPostProcessor") << "cannot find directory: " << subDir_ << " , skipping";
    return;
  }

  //loop over all triggers/samples//

  std::vector<std::string> subdirectories = dqm->getSubdirs();
  for(std::vector<std::string>::iterator dir = subdirectories.begin() ;dir!= subdirectories.end(); dir++ ){
    dqm->cd(*dir);


    // Do everything twice: once for mc-matched histos,   //
    // once for unmatched histos                          //

    std::vector<std::string> postfixes;
    std::string tmpstring=""; //unmatched histos
    postfixes.push_back(tmpstring);
    tmpstring="_MC_matched";
    postfixes.push_back(tmpstring);

    for(std::vector<std::string>::iterator postfix=postfixes.begin(); postfix!=postfixes.end();postfix++){
      
      // computer per-event efficiencies //
      
      std::string histoName="efficiency_by_step"+ *postfix;
      std::string baseName = "total_eff"+ *postfix;
      MonitorElement* total = dqm->book1D(histoName.c_str(),dqm->get(dqm->pwd() + "/" + baseName)->getTH1F());
      total->setTitle(histoName);
      
      //     std::vector<std::string> mes = dqm->getMEs();
      //     for(std::vector<std::string>::iterator me = mes.begin() ;me!= mes.end(); me++ )
      //       std::cout <<*me <<std::endl;
      //     std::cout <<std::endl;
      
      float value=0;
      float error=0;    
      //compute stepwise total efficiencies 
      for(int bin= total->getNbinsX()-2 ; bin > 1  ; bin--){
    value=0;
    error=0;
    if(total->getBinContent(bin-1) != 0){
      value = total->getBinContent(bin)/total->getBinContent(bin-1) ;
      error = sqrt(value*(1-value)/total->getBinContent(bin-1));
    }
    total->setBinContent(bin,value);
    total->setBinError(bin,error);
      }
      
      //set first bin to L1 efficiency
      if(total->getBinContent(total->getNbinsX()) !=0 ){
        value = total->getBinContent(1)/total->getBinContent(total->getNbinsX()) ;
    error = sqrt(value*(1-value)/total->getBinContent(total->getNbinsX()));
      }else{
      value=0;error=0;
      }
      total->setBinContent(1,value);
      total->setBinError(1,error);
      
      //total efficiency relative to gen
      if(total->getBinContent(total->getNbinsX()) !=0 ){
        value = dqm->get(dqm->pwd() + "/" + baseName)->getBinContent(total->getNbinsX()-2)/total->getBinContent(total->getNbinsX()) ;
    error = sqrt(value*(1-value)/total->getBinContent(total->getNbinsX()));
      }else{
    value=0;error=0;
      }
      total->setBinContent(total->getNbinsX(),value);
      total->setBinError(total->getNbinsX(),error);
      total->setBinLabel(total->getNbinsX(),"total efficiency rel. gen");
      
      //total efficiency relative to L1
      if(total->getBinContent(total->getNbinsX()) !=0 ){
        value = dqm->get(dqm->pwd() + "/" + baseName)->getBinContent(total->getNbinsX()-2)/dqm->get(dqm->pwd() + "/" + baseName)->getBinContent(1) ;
    error = sqrt(value*(1-value)/dqm->get(dqm->pwd() + "/" + baseName)->getBinContent(1));
      }else{
    value=0;error=0;
      }
      total->setBinContent(total->getNbinsX()-1,value);
      total->setBinError(total->getNbinsX()-1,error);
      total->setBinLabel(total->getNbinsX()-1,"total efficiency rel. L1");
      
      total->getTH1F()->SetMaximum(1.2);
      total->getTH1F()->SetMinimum(0);
      
      // compute per-object efficiencies       //
      MonitorElement *eff, *num, *denom, *genPlot, *effVsGen, *effL1VsGen;
      std::vector<std::string> varNames; 
      varNames.push_back("eta"); 
      varNames.push_back("phi"); 
      varNames.push_back("et");

      std::string filterName;
      std::string filterName2;
      std::string denomName;
      std::string numName;

      // Get the gen-level plots
      std::string genName;

      // Get the L1 over gen filter first
      filterName2= total->getTH1F()->GetXaxis()->GetBinLabel(1);
    
      //loop over variables (eta/phi/et)
      for(std::vector<std::string>::iterator var = varNames.begin(); var != varNames.end() ; var++){
    numName   = dqm->pwd() + "/" + filterName2 + *var + *postfix;
    genName   = dqm->pwd() + "/gen_" + *var ;
    num       = dqm->get(numName);
    genPlot   = dqm->get(genName);
        effL1VsGen = dqm->book1D("efficiency_"+filterName2+"_vs_"+*var +*postfix, dqm->get(numName)->getTH1F());

    // Check if histograms actually exist
    if(!num || !genPlot) break; 
    
    // Make sure we are able to book new element
    if (!dqm) break;

    // Create the efficiency plot
    effL1VsGen->setTitle("efficiency_"+filterName2+"_vs_"+*var + *postfix);
    effL1VsGen->getTH1F()->SetMaximum(1.2);
    effL1VsGen->getTH1F()->SetMinimum(0);
    effL1VsGen->getTH1F()->GetXaxis()->SetTitle(var->c_str());
    effL1VsGen->getTH1F()->Divide(num->getTH1F(),genPlot->getTH1F(),1,1,"b" );
      }
    
      // get the filter names from the bin-labels of the master-histogram
      for (int filter=1; filter < total->getNbinsX()-2; filter++) {
    filterName = total->getTH1F()->GetXaxis()->GetBinLabel(filter);
    filterName2= total->getTH1F()->GetXaxis()->GetBinLabel(filter+1);
    
    //loop over variables (eta/et)
    for(std::vector<std::string>::iterator var = varNames.begin(); var != varNames.end() ; var++){
      numName   = dqm->pwd() + "/" + filterName2 + *var + *postfix;
      denomName = dqm->pwd() + "/" + filterName  + *var + *postfix;
      num       = dqm->get(numName);
      denom     = dqm->get(denomName);

          // Check if histograms actually exist
      if(!num || !denom) break; 

      // Make sure we are able to book new element
          if (!dqm) break;

      // Is this the last filter? Book efficiency vs gen level
      std::string temp = *postfix;
          if (filter==total->getNbinsX()-3 && temp.find("matched")!=std::string::npos) {
            genName = dqm->pwd() + "/gen_" + *var;
            genPlot = dqm->get(genName);
            effVsGen = dqm->book1D("final_eff_vs_"+*var, dqm->get(genName)->getTH1F());
            if (!dqm) break;

            effVsGen->setTitle("Efficiency Compared to Gen vs "+*var);
        effVsGen->getTH1F()->SetMaximum(1.2); effVsGen->getTH1F()->SetMinimum(0.0);
        effVsGen->getTH1F()->GetXaxis()->SetTitle(var->c_str());
        effVsGen->getTH1F()->Divide(num->getTH1F(),genPlot->getTH1F(),1,1,"b" );
          }

      eff = dqm->book1D("efficiency_"+filterName2+"_vs_"+*var +*postfix, dqm->get(numName)->getTH1F());
      
          // Make sure we were able to book new element
          if (!dqm) break;

          // Create the efficiency plot
          /* num->getTH1F()->Sumw2();
      denom->getTH1F()->Sumw2(); 
          eff->getTH1F()->Sumw2(); */

      eff->setTitle("efficiency_"+filterName2+"_vs_"+*var + *postfix);
      eff->getTH1F()->SetMaximum(1.2);
      eff->getTH1F()->SetMinimum(0);
      eff->getTH1F()->GetXaxis()->SetTitle(var->c_str());
      eff->getTH1F()->Divide(num->getTH1F(),denom->getTH1F(),1,1,"b" );
    }
      }
    } 
    dqm->goUp();
  }
  
}
DEFINE_FWK_MODULE(EmDQMPostProcessor);