FSQDiJetAve.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    DQMOffline/FSQDiJetAve
// Class:      FSQDiJetAve
// 
//
// Original Author:  Tomasz Fruboes
//         Created:  Tue, 04 Nov 2014 11:36:27 GMT
//
//
// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"

#include "DQMOffline/Trigger/interface/FSQDiJetAve.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
#include "CommonTools/Utils/interface/StringObjectFunction.h"

#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include <DataFormats/TrackReco/interface/Track.h>
#include <DataFormats/EgammaCandidates/interface/Photon.h>
#include <DataFormats/MuonReco/interface/Muon.h>
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

#include <boost/algorithm/string.hpp>
#include "FWCore/Utilities/interface/EDGetToken.h"
#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"

using namespace edm;
using namespace std;

namespace FSQ {
//################################################################################################
//
// Base Handler class
//
//################################################################################################
class BaseHandler {
    public:
        BaseHandler();
        ~BaseHandler(){
        }
        BaseHandler(const edm::ParameterSet& iConfig,  triggerExpression::Data & eventCache):
              m_expression(triggerExpression::parse( iConfig.getParameter<std::string>("triggerSelection")))
         {
              // extract list of used paths
              std::vector<std::string> strs;
              std::string triggerSelection = iConfig.getParameter<std::string>("triggerSelection");
              boost::split(strs, triggerSelection, boost::is_any_of("\t ,`!@#$%^&*()~/\\"));
              for (size_t iToken = 0; iToken < strs.size();++iToken ){
                    if (strs.at(iToken).find("HLT_")==0){
                        m_usedPaths.insert(strs.at(iToken));
                    }
              }

              m_eventCache = &eventCache;
              std::string pathPartialName  = iConfig.getParameter<std::string>("partialPathName");
              m_dirname = iConfig.getUntrackedParameter("mainDQMDirname",std::string("HLT/FSQ/"))+pathPartialName + "/";
              m_pset = iConfig;

        };
        virtual void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
                             const HLTConfigProvider&  hltConfig,
                             const trigger::TriggerEvent& trgEvent,
                             const edm::TriggerResults & triggerResults, 
                             const edm::TriggerNames  & triggerNames,
                             float weight) = 0;
        virtual void book(DQMStore::IBooker & booker) = 0;
        virtual void getAndStoreTokens(edm::ConsumesCollector && iC) = 0;

        std::unique_ptr<triggerExpression::Evaluator>  m_expression;
        triggerExpression::Data * m_eventCache;
        std::string m_dirname;
        std::map<std::string,  MonitorElement*> m_histos;
        std::set<std::string> m_usedPaths;
        edm::ParameterSet m_pset;
};
//################################################################################################
//
// Handle objects saved into hlt event by hlt filters
//
//################################################################################################
enum SpecialFilters { None, BestVertexMatching, ApplyJEC };
template <class TInputCandidateType, class TOutputCandidateType, SpecialFilters filter = None>
class HandlerTemplate: public BaseHandler {
    private:
        std::string m_dqmhistolabel;
        std::string m_pathPartialName; //#("HLT_DiPFJetAve30_HFJEC_");
        std::string m_filterPartialName; //#("ForHFJECBase"); // Calo jet preFilter

        int m_combinedObjectDimension;

        StringCutObjectSelector<TInputCandidateType>  m_singleObjectSelection;
        StringCutObjectSelector<std::vector<TOutputCandidateType> >  m_combinedObjectSelection;
        StringObjectFunction<std::vector<TOutputCandidateType> >     m_combinedObjectSortFunction;
        std::map<std::string, std::shared_ptr<StringObjectFunction<std::vector<TOutputCandidateType> > > > m_plottersCombinedObject;
        std::map<std::string, std::shared_ptr<StringObjectFunction<TInputCandidateType > > > m_plottersSingleObject;
        static const int SingleObjectPlotter = 0;
        static const int CombinedObjectPlotter = 1;
        std::map<std::string, int > m_plotterType;
        std::vector< edm::ParameterSet > m_combinedObjectDrawables;
        std::vector< edm::ParameterSet > m_singleObjectDrawables; // for all single objects passing preselection
        bool m_isSetup;
        edm::InputTag m_input;
        std::map<std::string, edm::EDGetToken> m_tokens;

    public:
        HandlerTemplate(const edm::ParameterSet& iConfig, triggerExpression::Data & eventCache):
            BaseHandler(iConfig, eventCache),
            m_singleObjectSelection(iConfig.getParameter<std::string>("singleObjectsPreselection")),
            m_combinedObjectSelection(iConfig.getParameter<std::string>("combinedObjectSelection")),
            m_combinedObjectSortFunction(iConfig.getParameter<std::string>("combinedObjectSortCriteria"))
        {
             std::string type = iConfig.getParameter<std::string>("handlerType");
             if (type != "FromHLT") {
                m_input = iConfig.getParameter<edm::InputTag>("inputCol");
             }

             m_dqmhistolabel = iConfig.getParameter<std::string>("dqmhistolabel");
             m_filterPartialName = iConfig.getParameter<std::string>("partialFilterName"); // std::string find is used to match filter
             m_pathPartialName  = iConfig.getParameter<std::string>("partialPathName");
             m_combinedObjectDimension = iConfig.getParameter<int>("combinedObjectDimension");
             m_combinedObjectDrawables = iConfig.getParameter<  std::vector< edm::ParameterSet > >("combinedObjectDrawables");
             m_singleObjectDrawables = iConfig.getParameter<  std::vector< edm::ParameterSet > >("singleObjectDrawables");
             m_isSetup = false;
        }

        void book(DQMStore::IBooker & booker){
            if(!m_isSetup){
                booker.setCurrentFolder(m_dirname);
                m_isSetup = true;
                std::vector< std::vector< edm::ParameterSet > * > todo(2, (std::vector< edm::ParameterSet > * )0);
                todo[CombinedObjectPlotter]=&m_combinedObjectDrawables;
                todo[SingleObjectPlotter]=&m_singleObjectDrawables;
                for (size_t ti =0; ti<todo.size();++ti){
                    for (size_t i = 0; i < todo[ti]->size(); ++i){
                        std::string histoName = m_dqmhistolabel + "_" + todo[ti]->at(i).template getParameter<std::string>("name");
                        std::string expression =   todo[ti]->at(i).template getParameter<std::string>("expression");
                        int bins =    todo[ti]->at(i).template getParameter<int>("bins");
                        double rangeLow  =  todo[ti]->at(i).template getParameter<double>("min");
                        double rangeHigh =  todo[ti]->at(i).template getParameter<double>("max");
                        m_histos[histoName] =  booker.book1D(histoName, histoName, bins, rangeLow, rangeHigh);
                        m_plotterType[histoName] = ti;
                        if (ti == CombinedObjectPlotter){
                            StringObjectFunction<std::vector<TOutputCandidateType> > * func 
                                    = new StringObjectFunction<std::vector<TOutputCandidateType> >(expression);
                            m_plottersCombinedObject[histoName] =  std::shared_ptr<StringObjectFunction<std::vector<TOutputCandidateType> > >(func);
                        } else {
                            StringObjectFunction< TInputCandidateType>  * func 
                                    = new StringObjectFunction< TInputCandidateType> (expression);
                            m_plottersSingleObject[histoName] =  std::shared_ptr<StringObjectFunction<TInputCandidateType> > (func);
                        }
                    }   
                }
            }
        }
        void getAndStoreTokens(edm::ConsumesCollector && iC){
                edm::EDGetTokenT<std::vector<TInputCandidateType>  > tok =  iC.consumes<std::vector<TInputCandidateType> > (m_input);
                m_tokens[m_input.encode()] = edm::EDGetToken(tok);
        }

        //#############################################################################
        // Count objects. To avoid code duplication we do it in a separate template -
        //  - partial specialization not easy...:
        // http://stackoverflow.com/questions/21182729/specializing-single-method-in-a-big-template-class
        //#############################################################################
        template <class T>
        int count(const edm::Event& iEvent, InputTag &input, StringCutObjectSelector<T> & sel, float weight){
           int ret = 0;
           Handle<std::vector< T > > hIn;
           iEvent.getByToken(m_tokens[input.encode()], hIn);
           if(!hIn.isValid()) {
              edm::LogError("FSQDiJetAve") << "product not found: "<<  input.encode();
              return -1;  // return nonsense value
           }
           for (size_t i = 0; i<hIn->size(); ++i) {
                bool preselection = sel(hIn->at(i));
                if (preselection){
                    fillSingleObjectPlots(hIn->at(i), weight);
                    ret+=1;
                }
           }
           return ret;
        }

        // FIXME (?): code duplication
        void fillSingleObjectPlots(const TInputCandidateType & cand, float weight){
            std::map<std::string,  MonitorElement*>::iterator it, itE;
            it = m_histos.begin();
            itE = m_histos.end();
            for (;it!=itE;++it){
                if (m_plotterType[it->first]!=SingleObjectPlotter) continue;
                float val = (*m_plottersSingleObject[it->first])(cand);
                it->second->Fill(val, weight);
            }
        }
        // Notes:
        //  - this function (and specialized versions) are responsible for calling
        //     fillSingleObjectPlots for all single objects passing the single
        //     object preselection criteria
        //  - FIXME this function should take only event/ event setup (?)
        //  - FIXME responsibility to apply preselection should be elsewhere
        //          hard to fix, since we dont want to copy all objects due to
        //          performance reasons
        //  - note: implementation below working when in/out types are equal
        //          in other cases you must provide specialized version (see below)
        void getFilteredCands(TInputCandidateType *, std::vector<TOutputCandidateType> & cands,
                     const edm::Event& iEvent,
                     const edm::EventSetup& iSetup,
                     const HLTConfigProvider&  hltConfig,
                     const trigger::TriggerEvent& trgEvent,
                     float weight)
        {
               Handle<std::vector<TInputCandidateType> > hIn;
               iEvent.getByToken(m_tokens[m_input.encode()], hIn);

               if(!hIn.isValid()) {
                  edm::LogError("FSQDiJetAve") << "product not found: "<<  m_input.encode();
                  return;  
               }

               for (size_t i = 0; i<hIn->size(); ++i) {
                    bool preselection = m_singleObjectSelection(hIn->at(i));
                    if (preselection){
                        fillSingleObjectPlots(hIn->at(i), weight);
                        cands.push_back(hIn->at(i));
                    }
               }
        }

        std::vector<std::string > findPathAndFilter(const HLTConfigProvider&  hltConfig){
            std::vector<std::string> ret(2,"");
            std::string filterFullName = "";
            std::string pathFullName = "";
            std::vector<std::string> filtersForThisPath;
            //int pathIndex = -1;
            int numPathMatches = 0;
            int numFilterMatches = 0;
            for (size_t i = 0; i < hltConfig.size(); ++i) {
                if (hltConfig.triggerName(i).find(m_pathPartialName) == std::string::npos) continue;
                pathFullName = hltConfig.triggerName(i);
                //pathIndex = i;
                ++numPathMatches;
                std::vector<std::string > moduleLabels = hltConfig.moduleLabels(i);
                for (size_t iMod = 0; iMod <moduleLabels.size(); ++iMod){
                    if ("EDFilter" ==  hltConfig.moduleEDMType(moduleLabels.at(iMod))) {
                        filtersForThisPath.push_back(moduleLabels.at(iMod));
                        if ( moduleLabels.at(iMod).find(m_filterPartialName)!= std::string::npos  ){
                            filterFullName = moduleLabels.at(iMod);
                            ++numFilterMatches;
                        }
                    }
                }
            }

            // LogWarning or LogError?
            if (numPathMatches != 1) {
                  edm::LogError("FSQDiJetAve") << "Problem: found " << numPathMatches
                    << " paths matching " << m_pathPartialName << std::endl;
                  return ret;   
            }
            ret[0] = pathFullName;
            if (numFilterMatches != 1) {
                  edm::LogError("FSQDiJetAve") << "Problem: found " << numFilterMatches
                    << " filter matching " << m_filterPartialName
                    << " in path "<< m_pathPartialName << std::endl;
                  return ret;
            }
            ret[1] = filterFullName;
            return ret;
        }

        void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup,
                     const HLTConfigProvider&  hltConfig,
                     const trigger::TriggerEvent& trgEvent,
                     const edm::TriggerResults & triggerResults, 
                     const edm::TriggerNames  & triggerNames,
                     float weight)
        {
            size_t found = 0;
            for (size_t i = 0; i<triggerNames.size(); ++i){
                std::set<std::string>::iterator itUsedPaths = m_usedPaths.begin();
                for(; itUsedPaths != m_usedPaths.end(); ++itUsedPaths){ 
                    if (triggerNames.triggerName(i).find(*itUsedPaths)!= std::string::npos ){
                        ++found;
                        break;
                    }
                }

                if (found == m_usedPaths.size()) break;
            }
            if (found != m_usedPaths.size()){
                    edm::LogInfo("FSQDiJetAve") << "One of requested paths not found, skipping event";
                    return;
            }
            if (m_eventCache->configurationUpdated()) {
                m_expression->init(*m_eventCache);
            }
            if (not (*m_expression)(*m_eventCache)) return;

            /*
            std::vector<std::string> pathAndFilter = findPathAndFilter(hltConfig);

            std::string pathFullName = pathAndFilter[0];
            if (pathFullName == "") {
                return;
            }
            unsigned indexNum = triggerNames.triggerIndex(pathFullName);
            if(indexNum >= triggerNames.size()){
                  edm::LogError("FSQDiJetAve") << "Problem determining trigger index for " << pathFullName << " " << m_pathPartialName;
            }
            if (!triggerResults.accept(indexNum)) return;*/

            std::vector<TOutputCandidateType> cands;
            getFilteredCands((TInputCandidateType *)0, cands, iEvent, iSetup, hltConfig, trgEvent, weight);

            if (cands.size()==0) return;

            std::vector<TOutputCandidateType> bestCombinationFromCands = getBestCombination(cands);
            if (bestCombinationFromCands.size()==0) return;

            // plot 
            std::map<std::string,  MonitorElement*>::iterator it, itE;
            it = m_histos.begin();
            itE = m_histos.end();
            for (;it!=itE;++it){
                if (m_plotterType[it->first]!=CombinedObjectPlotter) continue;
                float val = (*m_plottersCombinedObject[it->first])(bestCombinationFromCands);
                it->second->Fill(val, weight);
            }
        }

        std::vector<TOutputCandidateType> getBestCombination(std::vector<TOutputCandidateType> & cands ){
            int columnSize = cands.size();
            std::vector<int> currentCombination(m_combinedObjectDimension, 0);
            std::vector<int> bestCombination(m_combinedObjectDimension, -1);

            int maxCombinations = 1;
            int cnt = 0;
            while (cnt < m_combinedObjectDimension){
                cnt += 1;
                maxCombinations *= columnSize;
            }

            cnt = 0;
            float bestCombinedCandVal = -1;
            while ( cnt < maxCombinations){
                cnt += 1;

                // 1. Check if current combination contains duplicates
                std::vector<int> currentCombinationCopy(currentCombination);
                std::vector<int>::iterator it;
                std::sort(currentCombinationCopy.begin(), currentCombinationCopy.end());
                it = std::unique(currentCombinationCopy.begin(), currentCombinationCopy.end());
                currentCombinationCopy.resize( std::distance(currentCombinationCopy.begin(),it) );
                bool duplicatesPresent = currentCombination.size() != currentCombinationCopy.size();
                // 2. If no duplicates found - 
                //          - check if current combination passes the cut
                //          - rank current combination
                if (!duplicatesPresent) { // no duplicates, we can consider this combined object
                    std::vector<TOutputCandidateType > currentCombinationFromCands;
                    for (int i = 0; i<m_combinedObjectDimension;++i){
                        currentCombinationFromCands.push_back( cands.at(currentCombination.at(i)));
                    }
                    bool isOK = m_combinedObjectSelection(currentCombinationFromCands);
                    if (isOK){
                        float curVal = m_combinedObjectSortFunction(currentCombinationFromCands);
                        // FIXME
                        if (curVal < 0) {
                            edm::LogError("FSQDiJetAve") << "Problem: ranking function returned negative value: " << curVal << std::endl;
                        } else if (curVal > bestCombinedCandVal){
                            //std::cout << curVal << " " << bestCombinedCandVal << std::endl;
                            bestCombinedCandVal = curVal;
                            bestCombination = currentCombination;
                        }
                    }
                }
                // 3. Prepare next combination to test
                //    note to future self: less error prone method with modulo
                currentCombination.at(m_combinedObjectDimension-1)+=1; // increase last number
                int carry = 0;
                for (int i = m_combinedObjectDimension-1; i>=0; --i){  // iterate over all numbers, check if we are out of range
                    currentCombination.at(i)+= carry;
                    carry = 0;
                    if (currentCombination.at(i)>=columnSize){
                        carry = 1;
                        currentCombination.at(i) = 0;
                    }
                }
            } // combinations loop ends

            std::vector<TOutputCandidateType > bestCombinationFromCands;
            if (bestCombination.size()!=0 && bestCombination.at(0)>=0){
                for (int i = 0; i<m_combinedObjectDimension;++i){
                          bestCombinationFromCands.push_back( cands.at(bestCombination.at(i)));
                }
            }
            return bestCombinationFromCands;
        }
};
//#############################################################################
// Read any object inheriting from reco::Candidate. Save p4
//
//  problem: for reco::Candidate there is no reflex dictionary, so selector
//  wont work
//#############################################################################
template<>
void HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector>::getAndStoreTokens(
                edm::ConsumesCollector && iC)
{
    edm::EDGetTokenT<View<reco::Candidate> > tok =  iC.consumes< View<reco::Candidate>  > (m_input);
    m_tokens[m_input.encode()] = edm::EDGetToken(tok);
}
template<>
void HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector>::getFilteredCands(
             reco::Candidate::LorentzVector *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<reco::Candidate::LorentzVector> & cands, // output collection
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent,
             float weight)
{  
   Handle<View<reco::Candidate> > hIn;
   iEvent.getByToken(m_tokens[m_input.encode()], hIn);
   if(!hIn.isValid()) {
      edm::LogError("FSQDiJetAve") << "product not found: "<<  m_input.encode();
      return;
   }
   for (size_t i = 0; i<hIn->size(); ++i) {
        bool preselection = m_singleObjectSelection(hIn->at(i).p4());
        if (preselection){
            fillSingleObjectPlots(hIn->at(i).p4(), weight);
            cands.push_back(hIn->at(i).p4());
        }
   }
}
//#############################################################################
//
// Count any object inheriting from reco::Track. Save into std::vector<int>
// note: this is similar to recoCand counter (code duplication is hard to 
//       avoid in this case)
//
//#############################################################################
template<>
void HandlerTemplate<reco::Track, int >::getFilteredCands(
             reco::Track *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<int > & cands, // output collection
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent,
             float weight)
{  
   cands.clear();
   cands.push_back(count<reco::Track>(iEvent, m_input, m_singleObjectSelection, weight) );
}
template<>
void HandlerTemplate<reco::GenParticle, int >::getFilteredCands(
             reco::GenParticle *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<int > & cands, // output collection
             const edm::Event& iEvent, const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig, const trigger::TriggerEvent& trgEvent, float weight)
{
   cands.clear();
   cands.push_back(count<reco::GenParticle>(iEvent, m_input, m_singleObjectSelection, weight) );
}
//#############################################################################
//
// Count any object inheriting from reco::Track that is not to distant from 
// selected vertex. Save into std::vector<int>
// note: this is similar to recoCand counter (code duplication is hard to 
//       avoid in this case)
//
//#############################################################################
template<>
void HandlerTemplate<reco::Track, int, BestVertexMatching>::getAndStoreTokens(
                edm::ConsumesCollector && iC)
{
    edm::EDGetTokenT<std::vector<reco::Track>  > tok =  iC.consumes<std::vector<reco::Track> > (m_input);
    m_tokens[m_input.encode()] = edm::EDGetToken(tok);

    edm::InputTag lVerticesTag = m_pset.getParameter<edm::InputTag>("vtxCollection");
    edm::EDGetTokenT<reco::VertexCollection > tok2 =  iC.consumes< reco::VertexCollection  > (lVerticesTag);
    m_tokens[lVerticesTag.encode()] = edm::EDGetToken(tok2);
}

template<>
void HandlerTemplate<reco::Track, int, BestVertexMatching>::getFilteredCands(
             reco::Track *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<int > & cands, // output collection
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent, float weight)
{  
    // this is not elegant, but should be thread safe
    static const edm::InputTag lVerticesTag = m_pset.getParameter<edm::InputTag>("vtxCollection");
    static const int lMinNDOF = m_pset.getParameter<int>("minNDOF"); //7
    static const double lMaxZ = m_pset.getParameter<double>("maxZ"); // 15
    static const double lMaxDZ = m_pset.getParameter<double>("maxDZ"); // 0.12
    static const double lMaxDZ2dzsigma = m_pset.getParameter<double>("maxDZ2dzsigma"); // 3
    static const double lMaxDXY = m_pset.getParameter<double>("maxDXY"); // 0.12
    static const double lMaxDXY2dxysigma = m_pset.getParameter<double>("maxDXY2dxysigma"); // 3

    cands.clear();
    cands.push_back(0);

    edm::Handle<reco::VertexCollection> vertices;
    iEvent.getByToken(m_tokens[lVerticesTag.encode()], vertices);

    //double bestvz=-999.9, bestvx=-999.9, bestvy=-999.9;

    double dxy, dz, dzsigma, dxysigma;
    math::XYZPoint vtxPoint(0.0,0.0,0.0);
    double vzErr =0.0, vxErr=0.0, vyErr=0.0;

    // take first vertex passing the criteria
    int bestVtx = -1;
    for (size_t i = 0; i < vertices->size(); ++i){
        if (vertices->at(i).ndof()<lMinNDOF) continue; 
        if (fabs(vertices->at(i).z())> lMaxZ) continue;

        vtxPoint=vertices->at(i).position();
        vzErr=vertices->at(i).zError();
        vxErr=vertices->at(i).xError();
        vyErr=vertices->at(i).yError();
        bestVtx = i;
        break;
    }
    if (bestVtx < 0) return;
    // const reco::Vertex & vtx = vertices->at(bestVtx);

   Handle<std::vector<reco::Track > > hIn;
   iEvent.getByToken(m_tokens[m_input.encode()], hIn);
   if(!hIn.isValid()) {
      edm::LogError("FSQDiJetAve") << "product not found: "<<  m_input.encode();
      return;
   }

   for (size_t i = 0; i<hIn->size(); ++i) {
        if (!m_singleObjectSelection(hIn->at(i))) continue;
        dxy=0.0, dz=0.0, dxysigma=0.0, dzsigma=0.0;
        dxy = -1.*hIn->at(i).dxy(vtxPoint);
        dz = hIn->at(i).dz(vtxPoint);
        dxysigma = sqrt(hIn->at(i).dxyError()*hIn->at(i).dxyError()+vxErr*vyErr);
        dzsigma = sqrt(hIn->at(i).dzError()*hIn->at(i).dzError()+vzErr*vzErr);
        
        if(fabs(dz)>lMaxDZ)continue; // TODO...
        if(fabs(dz/dzsigma)>lMaxDZ2dzsigma)continue;
        if(fabs(dxy)>lMaxDXY)continue;
        if(fabs(dxy/dxysigma)>lMaxDXY2dxysigma)continue;
        
        cands.at(0)+=1;
    }//loop over tracks
}
//#############################################################################
//
// Apply JEC to PFJets
//
//#############################################################################
template<>
void HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC>::getAndStoreTokens(
                edm::ConsumesCollector && iC)
{
    edm::EDGetTokenT<std::vector<reco::PFJet>  > tok =  iC.consumes<std::vector<reco::PFJet> > (m_input);
    m_tokens[m_input.encode()] = edm::EDGetToken(tok);

    edm::InputTag jetCorTag = m_pset.getParameter<edm::InputTag>("PFJetCorLabel");
    edm::EDGetTokenT<reco::JetCorrector> jetcortoken =  iC.consumes<reco::JetCorrector>(jetCorTag);
    m_tokens[jetCorTag.encode()] = edm::EDGetToken(jetcortoken);

}

template<>
void HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC>::getFilteredCands(
             reco::PFJet *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<reco::PFJet > & cands, // output collection
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent, float weight)
{  
    cands.clear();
    static const edm::InputTag jetCorTag = m_pset.getParameter<edm::InputTag>("PFJetCorLabel");
    edm::Handle<reco::JetCorrector> pfcorrector;
    iEvent.getByToken(m_tokens[jetCorTag.encode()], pfcorrector);
 
    Handle<std::vector<reco::PFJet> > hIn;
    iEvent.getByToken(m_tokens[m_input.encode()], hIn);

    if(!hIn.isValid()) {
      edm::LogError("FSQDiJetAve") << "product not found: "<<  m_input.encode();
      return;  
    }

    for (size_t i = 0; i<hIn->size(); ++i) {
         double scale = pfcorrector->correction(hIn->at(i));
         reco::PFJet newPFJet(scale*hIn->at(i).p4(), hIn->at(i).vertex(), 
                              hIn->at(i).getSpecific(),  hIn->at(i).getJetConstituents());
         
         bool preselection = m_singleObjectSelection(newPFJet);
         if (preselection){
             fillSingleObjectPlots(newPFJet, weight);
             cands.push_back(newPFJet);
         }
    }
}
//#############################################################################
//
// Count any object inheriting from reco::Candidate. Save into std::vector<int>
//  same problem as for reco::Candidate handler ()
//
//#############################################################################
template<>
void HandlerTemplate<reco::Candidate::LorentzVector, int>::getAndStoreTokens(
                edm::ConsumesCollector && iC)
{
    edm::EDGetTokenT<View<reco::Candidate> > tok =  iC.consumes< View<reco::Candidate>  > (m_input);
    m_tokens[m_input.encode()] = edm::EDGetToken(tok);
}
template<>
void HandlerTemplate<reco::Candidate::LorentzVector, int >::getFilteredCands(
             reco::Candidate::LorentzVector *, // pass a dummy pointer, makes possible to select correct getFilteredCands
             std::vector<int > & cands, // output collection
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent,
             float weight)
{  
   cands.clear();
   cands.push_back(0);

   Handle<View<reco::Candidate> > hIn;
   iEvent.getByToken(m_tokens[m_input.encode()], hIn);
   if(!hIn.isValid()) {
      edm::LogError("FSQDiJetAve") << "product not found: "<<  m_input.encode();
      return;
   }
   for (size_t i = 0; i<hIn->size(); ++i) {
        bool preselection = m_singleObjectSelection(hIn->at(i).p4());
        if (preselection){
            fillSingleObjectPlots(hIn->at(i).p4(), weight);
            cands.at(0)+=1;
        }
   }
}
//#############################################################################
//
// Read and save trigger::TriggerObject from triggerEvent
//
//#############################################################################
template<>
void HandlerTemplate<trigger::TriggerObject, trigger::TriggerObject>::getFilteredCands(
             trigger::TriggerObject *, 
             std::vector<trigger::TriggerObject> &cands, 
             const edm::Event& iEvent,  
             const edm::EventSetup& iSetup,
             const HLTConfigProvider&  hltConfig,
             const trigger::TriggerEvent& trgEvent,
             float weight)
{
    // 1. Find matching path. Inside matchin path find matching filter
    std::string filterFullName = findPathAndFilter(hltConfig)[1];
    if (filterFullName == "") {
        return;
    }

    // 2. Fetch HLT objects saved by selected filter. Save those fullfilling preselection
    //      objects are saved in cands variable
    std::string process = trgEvent.usedProcessName(); // broken?
    edm::InputTag hltTag(filterFullName ,"", process);
    
    const int hltIndex = trgEvent.filterIndex(hltTag);
    if ( hltIndex >= trgEvent.sizeFilters() ) {
      edm::LogInfo("FSQDiJetAve") << "Cannot determine hlt index for |" << filterFullName << "|" << process;
      return;
    }

    const trigger::TriggerObjectCollection & toc(trgEvent.getObjects());
    const trigger::Keys & khlt = trgEvent.filterKeys(hltIndex);

    trigger::Keys::const_iterator kj = khlt.begin();

    for(;kj != khlt.end(); ++kj){
        bool preselection = m_singleObjectSelection(toc[*kj]);
        if (preselection){
            fillSingleObjectPlots(toc[*kj], weight);
            cands.push_back( toc[*kj]);
        }
    }

}

typedef HandlerTemplate<trigger::TriggerObject, trigger::TriggerObject> HLTHandler;
typedef HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector> RecoCandidateHandler;// in fact reco::Candidate, reco::Candidate::LorentzVector
typedef HandlerTemplate<reco::PFJet, reco::PFJet> RecoPFJetHandler;
typedef HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC> RecoPFJetWithJECHandler;
typedef HandlerTemplate<reco::Track, reco::Track> RecoTrackHandler;
typedef HandlerTemplate<reco::Photon, reco::Photon> RecoPhotonHandler;
typedef HandlerTemplate<reco::Muon, reco::Muon> RecoMuonHandler;
typedef HandlerTemplate<reco::GenParticle, reco::GenParticle > RecoGenParticleHandler;
typedef HandlerTemplate<reco::Candidate::LorentzVector, int > RecoCandidateCounter;
typedef HandlerTemplate<reco::Track, int > RecoTrackCounter;
typedef HandlerTemplate<reco::Track, int, BestVertexMatching> RecoTrackCounterWithVertexConstraint;
typedef HandlerTemplate<reco::GenParticle, int > RecoGenParticleCounter;
}
//################################################################################################
//
// Plugin functions
//
//################################################################################################
FSQDiJetAve::FSQDiJetAve(const edm::ParameterSet& iConfig):
  m_eventCache(iConfig.getParameterSet("triggerConfiguration") , consumesCollector())
{
  m_useGenWeight = iConfig.getParameter<bool>("useGenWeight");
  if (m_useGenWeight) {
       m_genEvInfoToken = consumes <GenEventInfoProduct> (edm::InputTag("generator"));
  }

  triggerSummaryLabel_ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
  triggerResultsLabel_ = iConfig.getParameter<edm::InputTag>("triggerResultsLabel");
  triggerSummaryToken  = consumes <trigger::TriggerEvent> (triggerSummaryLabel_);
  triggerResultsToken  = consumes <edm::TriggerResults>   (triggerResultsLabel_);

  triggerSummaryFUToken= consumes <trigger::TriggerEvent> (edm::InputTag(triggerSummaryLabel_.label(),triggerSummaryLabel_.instance(),std::string("FU")));
  triggerResultsFUToken= consumes <edm::TriggerResults>   (edm::InputTag(triggerResultsLabel_.label(),triggerResultsLabel_.instance(),std::string("FU")));

  std::vector< edm::ParameterSet > todo  = iConfig.getParameter<  std::vector< edm::ParameterSet > >("todo");
  for (size_t i = 0; i < todo.size(); ++i) {
        edm::ParameterSet pset = todo.at(i);
        std::string type = pset.getParameter<std::string>("handlerType");
        if (type == "FromHLT") {
            m_handlers.push_back(std::shared_ptr<FSQ::HLTHandler>(new FSQ::HLTHandler(pset, m_eventCache)));
        }
        else if (type == "RecoCandidateCounter") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoCandidateCounter>(new FSQ::RecoCandidateCounter(pset, m_eventCache)));
        }
        else if (type == "RecoTrackCounter") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoTrackCounter>(new FSQ::RecoTrackCounter(pset, m_eventCache)));
        }
        else if (type == "RecoTrackCounterWithVertexConstraint") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoTrackCounterWithVertexConstraint>
                    (new FSQ::RecoTrackCounterWithVertexConstraint(pset, m_eventCache)));
        }
        else if (type == "FromRecoCandidate") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoCandidateHandler>(new FSQ::RecoCandidateHandler(pset, m_eventCache)));
        }
        else if (type == "RecoPFJet") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoPFJetHandler>(new FSQ::RecoPFJetHandler(pset, m_eventCache)));
        }
        else if (type == "RecoPFJetWithJEC") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoPFJetWithJECHandler>(new FSQ::RecoPFJetWithJECHandler(pset, m_eventCache)));
        }
        else if (type == "RecoTrack") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoTrackHandler>(new FSQ::RecoTrackHandler(pset, m_eventCache)));
        } 
        else if (type == "RecoPhoton") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoPhotonHandler>(new FSQ::RecoPhotonHandler(pset, m_eventCache)));
        } 
        else if (type == "RecoMuon") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoMuonHandler>(new FSQ::RecoMuonHandler(pset, m_eventCache)));
        } 
        else if (type == "RecoGenParticleCounter") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoGenParticleCounter>(new FSQ::RecoGenParticleCounter(pset, m_eventCache)));
        }
        else if (type == "RecoGenParticleHandler") {
            m_handlers.push_back(std::shared_ptr<FSQ::RecoGenParticleHandler>(new FSQ::RecoGenParticleHandler(pset, m_eventCache)));
        } 
        else {
            throw cms::Exception("FSQ DQM handler not know: "+ type);
        }
  }
  for (size_t i = 0; i < m_handlers.size(); ++i) {
        m_handlers.at(i)->getAndStoreTokens(consumesCollector());
  }

}

FSQDiJetAve::~FSQDiJetAve()
{}

void
FSQDiJetAve::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;
   if (not m_eventCache.setEvent(iEvent, iSetup)){
      edm::LogError("FSQDiJetAve") << "Could not setup the filter";
   }

  //---------- triggerResults ----------
  iEvent.getByToken(triggerResultsToken, m_triggerResults);
  if(!m_triggerResults.isValid()) {
    iEvent.getByToken(triggerResultsFUToken, m_triggerResults);
    if(!m_triggerResults.isValid()) {
      edm::LogError("FSQDiJetAve") << "TriggerResults not valid, skippng event";
      return;
    }
  }
  
  //---------- triggerResults ----------
  if(m_triggerResults.isValid()) {
    m_triggerNames = iEvent.triggerNames(*m_triggerResults);
  } else {
    edm::LogError("FSQDiJetAve") << "TriggerResults not found";
    return;
  } 
  
  //---------- triggerSummary ----------
  iEvent.getByToken(triggerSummaryToken, m_trgEvent);
  if(!m_trgEvent.isValid()) {
    iEvent.getByToken(triggerSummaryFUToken, m_trgEvent);
    if(!m_trgEvent.isValid()) {
      edm::LogInfo("FSQDiJetAve") << "TriggerEvent not found, ";
      return;
    }
  } 

  float weight = 1.;  
  if (m_useGenWeight){
    edm::Handle<GenEventInfoProduct> hGW;
    iEvent.getByToken(m_genEvInfoToken, hGW);
    weight = hGW->weight();
  }

  for (size_t i = 0; i < m_handlers.size(); ++i) {
        m_handlers.at(i)->analyze(iEvent, iSetup, m_hltConfig, *m_trgEvent.product(), *m_triggerResults.product(), m_triggerNames, weight);
  }

}
// ------------ method called when starting to processes a run  ------------
//*
void 
FSQDiJetAve::dqmBeginRun(edm::Run const& run, edm::EventSetup const& c)
{
    bool changed(true);
    std::string processName = triggerResultsLabel_.process();
    if (m_hltConfig.init(run, c, processName, changed)) {
        LogDebug("FSQDiJetAve") << "HLTConfigProvider failed to initialize.";
    }

}
void FSQDiJetAve::bookHistograms(DQMStore::IBooker & booker, edm::Run const & run, edm::EventSetup const & c){
    for (size_t i = 0; i < m_handlers.size(); ++i) {
        m_handlers.at(i)->book(booker);
    }
}
//*/
// ------------ method called when ending the processing of a run  ------------
/*
void 
FSQDiJetAve::endRun(edm::Run const&, edm::EventSetup const&)
{
}
*/

// ------------ method called when starting to processes a luminosity block  ------------
/*
void 
FSQDiJetAve::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{
}
*/

// ------------ method called when ending the processing of a luminosity block  ------------
/*
void 
FSQDiJetAve::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{}
// */

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void
FSQDiJetAve::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}