HLTHiggsPlotter.cc

Go to the documentation of this file.

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Candidate/interface/CandMatchMap.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"

#include "DataFormats/HLTReco/interface/TriggerEvent.h"

#include "HLTriggerOffline/Higgs/interface/HLTHiggsPlotter.h"
#include "HLTriggerOffline/Higgs/interface/HLTHiggsSubAnalysis.h"
#include "HLTriggerOffline/Higgs/src/EVTColContainer.cc"

#include "TPRegexp.h"


#include<set>
#include<cctype>

HLTHiggsPlotter::HLTHiggsPlotter(const edm::ParameterSet & pset,
                                 const std::string & hltPath,
                                 const std::vector<unsigned int> & objectsType, 
                                 const unsigned int & NptPlots,
                                 const std::vector<double> & NminOneCuts) :
    _hltPath(hltPath),
    _hltProcessName(pset.getParameter<std::string>("hltProcessName")),
    _objectsType(std::set<unsigned int>(objectsType.begin(),objectsType.end())),
    _nObjects(objectsType.size()),
    _parametersEta(pset.getParameter<std::vector<double> >("parametersEta")),
    _parametersPhi(pset.getParameter<std::vector<double> >("parametersPhi")),
    _parametersTurnOn(pset.getParameter<std::vector<double> >("parametersTurnOn")),
    _NptPlots(NptPlots),
    _NminOneCuts(NminOneCuts)
{
    for(std::set<unsigned int>::iterator it = _objectsType.begin();
        it != _objectsType.end(); ++it)
    {
        // Some parameters extracted from the .py
        std::string objStr = EVTColContainer::getTypeString( *it );
        _cutMinPt[*it] = pset.getParameter<double>( std::string(objStr+"_cutMinPt").c_str() );
        _cutMaxEta[*it] = pset.getParameter<double>( std::string(objStr+"_cutMaxEta").c_str() );
    }
}

HLTHiggsPlotter::~HLTHiggsPlotter()
{
}


void HLTHiggsPlotter::beginJob() 
{
}



void HLTHiggsPlotter::beginRun(const edm::Run & iRun,
                            const edm::EventSetup & iSetup)
{

}

void HLTHiggsPlotter::bookHistograms(DQMStore::IBooker &ibooker, const bool & useNminOneCuts)
{
    for (std::set<unsigned int>::iterator it = _objectsType.begin();
        it != _objectsType.end(); ++it)
    {
        std::vector<std::string> sources(2);
        sources[0] = "gen";
        sources[1] = "rec";
        TString maxPt;
        
        const std::string objTypeStr = EVTColContainer::getTypeString(*it);
        
        for (size_t i = 0; i < sources.size(); i++)
        {
            std::string source = sources[i];
            
            if( useNminOneCuts && *it == EVTColContainer::PFJET) {
                if( source == "gen" ) continue;
                else {
                    // N-1 jet plots (dEtaqq, mqq, dPhibb, CSV1, maxCSV_jets, maxCSV_E, PFMET, pt1, pt2, pt3, pt4)
                    if( _NminOneCuts[0] ) bookHist(source, objTypeStr, "dEtaqq", ibooker);
                    if( _NminOneCuts[1] ) bookHist(source, objTypeStr, "mqq", ibooker);
                    if( _NminOneCuts[2] ) bookHist(source, objTypeStr, "dPhibb", ibooker);
                    if( _NminOneCuts[3] ) {
                        if ( _NminOneCuts[6] ) bookHist(source, objTypeStr, "maxCSV", ibooker);
                        else bookHist(source, objTypeStr, "CSV1", ibooker);
                    }
                    if( _NminOneCuts[4] ) bookHist(source, objTypeStr, "CSV2", ibooker);
                    if( _NminOneCuts[5] ) bookHist(source, objTypeStr, "CSV3", ibooker);
                }
            }
            
            bookHist(source, objTypeStr, "Eta", ibooker);
            bookHist(source, objTypeStr, "Phi", ibooker);
            for( unsigned int i=0; i < _NptPlots; i++ )
            {
                maxPt = "MaxPt";
                maxPt += i+1;
                bookHist(source, objTypeStr, maxPt.Data(), ibooker);
            }
        }
    }  
}

void HLTHiggsPlotter::analyze(const bool & isPassTrigger,
                            const std::string & source,
                            const std::vector<MatchStruct> & matches,
                            const unsigned int & minCandidates
                            )
{
    if ( !isPassTrigger )
    {
        return;
    }
    std::map<unsigned int,int> countobjects;
    // Initializing the count of the used object
    for(std::set<unsigned int>::iterator co = _objectsType.begin();
        co != _objectsType.end(); ++co)
    {
        countobjects[*co] = 0;
    }
        
    int counttotal = 0;
    const int totalobjectssize2 = _NptPlots*countobjects.size();
    // Fill the histos if pass the trigger (just the two with higher pt)
    for (size_t j = 0; j < matches.size(); ++j)
    {
        // Is this object owned by this trigger? If not we are not interested...
        if ( _objectsType.find( matches[j].objType) == _objectsType.end() )
        {
            continue;
        }

        const unsigned int objType = matches[j].objType;
        const std::string objTypeStr = EVTColContainer::getTypeString(matches[j].objType);
            
        float pt  = matches[j].pt;
        float eta = matches[j].eta;
        float phi = matches[j].phi;
        
        TString maxPt;
        if( (unsigned)countobjects[objType] < _NptPlots )
        {
            maxPt = "MaxPt";
            maxPt += (countobjects[objType]+1);
            this->fillHist(isPassTrigger,source,objTypeStr,maxPt.Data(),pt);
            // Filled the high pt ...
            ++(countobjects[objType]);
            ++counttotal;
        }
        else {
            if( (unsigned)countobjects[objType] < minCandidates ) { // To get correct results for HZZ
                ++(countobjects[objType]);
                ++counttotal;
            }
            else continue; //   Otherwise too many entries in Eta and Phi distributions
        }
        
        this->fillHist(isPassTrigger,source,objTypeStr,"Eta",eta);
        this->fillHist(isPassTrigger,source,objTypeStr,"Phi",phi);

        if ( counttotal == totalobjectssize2 ) 
        {
            break;
        }               
    }
}

void HLTHiggsPlotter::analyze(const bool & isPassTrigger, const std::string & source, const std::vector<MatchStruct> & matches,
                std::map<std::string,bool> & nMinOne, const float & dEtaqq, const float & mqq, const float & dPhibb, const float & CSV1, const float & CSV2, const float & CSV3, const bool & passAllCuts)
{
    if ( !isPassTrigger )
    {
        return;
    }
    std::map<unsigned int,int> countobjects;
    // Initializing the count of the used object
    for(std::set<unsigned int>::iterator co = _objectsType.begin();
        co != _objectsType.end(); ++co)
    {
        if( !(*co == EVTColContainer::PFJET && source == "gen") ) // genJets are not there
            countobjects[*co] = 0;
    }
        
    int counttotal = 0;
    const int totalobjectssize2 = _NptPlots*countobjects.size();
    // Fill the histos if pass the trigger (just the two with higher pt)
    for (size_t j = 0; j < matches.size(); ++j)
    { 
        // Is this object owned by this trigger? If not we are not interested...
        if ( _objectsType.find( matches[j].objType) == _objectsType.end() )
        {
            continue;
        }

        const unsigned int objType = matches[j].objType;
        const std::string objTypeStr = EVTColContainer::getTypeString(matches[j].objType);
            
        float pt  = matches[j].pt;
        float eta = matches[j].eta;
        float phi = matches[j].phi;
        
        // PFMET N-1 cut
        if( objType == EVTColContainer::PFMET && _NminOneCuts[8] && ! nMinOne["PFMET"] ) continue;
        
        TString maxPt;
        if( (unsigned)(countobjects)[objType] < _NptPlots )
        {
            maxPt = "MaxPt";
            maxPt += (countobjects[objType]+1);
            if( objType != EVTColContainer::PFJET || nMinOne[maxPt.Data()] ) {
                this->fillHist(isPassTrigger,source,objTypeStr,maxPt.Data(),pt);
            }
            ++(countobjects[objType]);
            ++counttotal;
        } 
        else continue; // if not needed (minCandidates == _NptPlots if _useNminOneCuts 
        if( objType != EVTColContainer::PFJET || passAllCuts ) {
            this->fillHist(isPassTrigger,source,objTypeStr,"Eta",eta);
            this->fillHist(isPassTrigger,source,objTypeStr,"Phi",phi);
        }
        
        if ( counttotal == totalobjectssize2 ) 
        {
            break;
        }
    }
    if( source == "rec" && _objectsType.find(EVTColContainer::PFJET) != _objectsType.end()) {
        if( _NminOneCuts[0] && nMinOne["dEtaqq"] ) {
            this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),"dEtaqq",dEtaqq);
        }
        if( _NminOneCuts[1] && nMinOne["mqq"] ) {
            this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),"mqq",mqq);
        }
        if( _NminOneCuts[2] && nMinOne["dPhibb"] ) {
            this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),"dPhibb",dPhibb);
        }
        if( _NminOneCuts[3] ) {
            std::string nameCSVplot = "CSV1";
            if ( _NminOneCuts[6] ) nameCSVplot = "maxCSV";
            if ( nMinOne[nameCSVplot] ) this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),nameCSVplot,CSV1);
        }
        if( _NminOneCuts[4] && nMinOne["CSV2"] ) {
            this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),"CSV2",CSV2);
        }
        if( _NminOneCuts[5] && nMinOne["CSV3"] ) {
            this->fillHist(isPassTrigger,source,EVTColContainer::getTypeString(EVTColContainer::PFJET),"CSV3",CSV3);
        }
    }
}

void HLTHiggsPlotter::bookHist(const std::string & source, 
                            const std::string & objType,
                            const std::string & variable,
                            DQMStore::IBooker & ibooker)
{
    std::string sourceUpper = source; 
    sourceUpper[0] = std::toupper(sourceUpper[0]);
    std::string name = source + objType + variable + "_" + _hltPath;
    TH1F * h = nullptr;

    if (variable.find("MaxPt") != std::string::npos) 
    {
        std::string desc;
        if (variable == "MaxPt1") desc = "Leading";
        else if (variable == "MaxPt2") desc = "Next-to-Leading";
        else desc = variable.substr(5,6) + "th Leading";
        std::string title = "pT of " + desc + " " + sourceUpper + " " + objType + " "
            "where event pass the "+ _hltPath;
        const size_t nBinsStandard = _parametersTurnOn.size() - 1;
        size_t nBins = nBinsStandard;
        float * edges = new float[nBinsStandard + 1];
        for(size_t i = 0; i < nBinsStandard + 1; i++)
        {
            edges[i] = _parametersTurnOn[i];
        }
        
        std::string jetObj = EVTColContainer::getTypeString(EVTColContainer::PFJET);
        if( objType == jetObj ) {
            const size_t nBinsJets = 25;
            nBins = nBinsJets;
            delete [] edges;
            edges = new float[nBinsJets+1];
            for(size_t i = 0; i < nBinsJets + 1; i++)
            {
                edges[i] = i*10;
            }
        }
        if( objType == EVTColContainer::getTypeString(EVTColContainer::PFMET) ) {
            const size_t nBinsJets = 30;
            nBins = nBinsJets;
            delete [] edges;
            edges = new float[nBinsJets+1];
            for(size_t i = 0; i < nBinsJets + 1; i++)
            {
                edges[i] = i*10;
            }
        }
        h = new TH1F(name.c_str(), title.c_str(), nBins, edges);
        delete [] edges;
    }
    else 
    {
        if( variable == "dEtaqq" ){
            std::string title  = "#Delta #eta_{qq} of " + sourceUpper + " " + objType;
            int    nBins = 20;
            double min   = 0;
            double max   = 4.8;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        }
        else if ( variable == "mqq" ){
            std::string title  = "m_{qq} of " + sourceUpper + " " + objType;
            int nBins = 20;
            double min   = 0;
            double max   = 1000;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        } 
        else if ( variable == "dPhibb" ){
            std::string title  = "#Delta #phi_{bb} of " + sourceUpper + " " + objType;
            int    nBins = 10;
            double min   = 0;
            double max   = 3.1416;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        }
        else if ( variable == "CSV1" ){
            std::string title  = "CSV1 of " + sourceUpper + " " + objType;
            int    nBins = 20;
            double min   = 0;
            double max   = 1;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        } 
        else if ( variable == "CSV2" ){
            std::string title  = "CSV2 of " + sourceUpper + " " + objType;
            int    nBins = 20;
            double min   = 0;
            double max   = 1;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        } 
        else if ( variable == "CSV3" ){
            std::string title  = "CSV3 of " + sourceUpper + " " + objType;
            int    nBins = 20;
            double min   = 0;
            double max   = 1;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        } 
        else if ( variable == "maxCSV" ){
            std::string title  = "max CSV of " + sourceUpper + " " + objType;
            int    nBins = 20;
            double min   = 0;
            double max   = 1;
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        } 
        else
        {
            std::string symbol = (variable == "Eta") ? "#eta" : "#phi";
            std::string title  = symbol + " of " + sourceUpper + " " + objType + " "+
            "where event pass the "+ _hltPath;
            std::vector<double> params = (variable == "Eta") ? _parametersEta : _parametersPhi;

            int    nBins = (int)params[0];
            double min   = params[1];
            double max   = params[2];
            h = new TH1F(name.c_str(), title.c_str(), nBins, min, max);
        }
    }
    h->Sumw2();
    _elements[name] = ibooker.book1D(name.c_str(), h);
    delete h;
}

void HLTHiggsPlotter::fillHist(const bool & passTrigger,
                               const std::string & source, 
                               const std::string & objType,
                               const std::string & variable,
                               const float & value )
{
    std::string sourceUpper = source; 
    sourceUpper[0] = toupper(sourceUpper[0]);
    std::string name = source + objType + variable + "_" + _hltPath;

    _elements[name]->Fill(value);
}