Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:33

0001 // -*- C++ -*-
0002 //
0003 // Package:    DQMOffline/FSQDiJetAve
0004 // Class:      FSQDiJetAve
0005 //
0006 /**\class FSQDiJetAve FSQDiJetAve.cc DQMOffline/FSQDiJetAve/plugins/FSQDiJetAve.cc
0007 
0008  Description: DQM source for FSQ triggers
0009 
0010  Implementation:
0011 */
0012 //
0013 // Original Author:  Tomasz Fruboes
0014 //         Created:  Tue, 04 Nov 2014 11:36:27 GMT
0015 //
0016 //
0017 // system include files
0018 #include <memory>
0019 
0020 // user include files
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 #include "FWCore/Framework/interface/ConsumesCollector.h"
0029 
0030 #include "DQMOffline/Trigger/interface/FSQDiJetAve.h"
0031 
0032 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0033 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0034 
0035 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0036 
0037 #include "DataFormats/Candidate/interface/Candidate.h"
0038 #include "DataFormats/JetReco/interface/PFJet.h"
0039 #include <DataFormats/TrackReco/interface/Track.h>
0040 #include <DataFormats/EgammaCandidates/interface/Photon.h>
0041 #include <DataFormats/MuonReco/interface/Muon.h>
0042 #include "DataFormats/VertexReco/interface/Vertex.h"
0043 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0044 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0045 
0046 #include <boost/algorithm/string.hpp>
0047 #include "FWCore/Utilities/interface/EDGetToken.h"
0048 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0049 
0050 using namespace edm;
0051 using namespace std;
0052 
0053 namespace FSQ {
0054   //################################################################################################
0055   //
0056   // Base Handler class
0057   //
0058   //################################################################################################
0059   class BaseHandler {
0060   public:
0061     typedef dqm::legacy::MonitorElement MonitorElement;
0062     typedef dqm::legacy::DQMStore DQMStore;
0063 
0064     BaseHandler();
0065     virtual ~BaseHandler() = default;
0066     BaseHandler(const edm::ParameterSet& iConfig, triggerExpression::Data& eventCache)
0067         : m_expression(triggerExpression::parse(iConfig.getParameter<std::string>("triggerSelection"))) {
0068       // extract list of used paths
0069       std::vector<std::string> strs;
0070       std::string triggerSelection = iConfig.getParameter<std::string>("triggerSelection");
0071       boost::split(strs, triggerSelection, boost::is_any_of("\t ,`!@#$%^&*()~/\\"));
0072       for (auto& str : strs) {
0073         if (str.find("HLT_") == 0) {
0074           m_usedPaths.insert(str);
0075         }
0076       }
0077 
0078       m_eventCache = &eventCache;
0079       std::string pathPartialName = iConfig.getParameter<std::string>("partialPathName");
0080       m_dirname = iConfig.getUntrackedParameter("mainDQMDirname", std::string("HLT/FSQ/")) + pathPartialName + "/";
0081       m_pset = iConfig;
0082     };
0083     virtual void analyze(const edm::Event& iEvent,
0084                          const edm::EventSetup& iSetup,
0085                          const HLTConfigProvider& hltConfig,
0086                          const trigger::TriggerEvent& trgEvent,
0087                          const edm::TriggerResults& triggerResults,
0088                          const edm::TriggerNames& triggerNames,
0089                          float weight) = 0;
0090     virtual void book(DQMStore::IBooker& booker) = 0;
0091     virtual void getAndStoreTokens(edm::ConsumesCollector&& iC) = 0;
0092 
0093     std::unique_ptr<triggerExpression::Evaluator> m_expression;
0094     triggerExpression::Data* m_eventCache;
0095     std::string m_dirname;
0096     std::map<std::string, MonitorElement*> m_histos;
0097     std::set<std::string> m_usedPaths;
0098     edm::ParameterSet m_pset;
0099   };
0100   //################################################################################################
0101   //
0102   // Handle objects saved into hlt event by hlt filters
0103   //
0104   //################################################################################################
0105   enum SpecialFilters { None, BestVertexMatching, ApplyJEC };
0106   template <class TInputCandidateType, class TOutputCandidateType, SpecialFilters filter = None>
0107   class HandlerTemplate : public BaseHandler {
0108   private:
0109     std::string m_dqmhistolabel;
0110     std::string m_pathPartialName;    //#("HLT_DiPFJetAve30_HFJEC_");
0111     std::string m_filterPartialName;  //#("ForHFJECBase"); // Calo jet preFilter
0112 
0113     int m_combinedObjectDimension;
0114 
0115     StringCutObjectSelector<TInputCandidateType> m_singleObjectSelection;
0116     StringCutObjectSelector<std::vector<TOutputCandidateType> > m_combinedObjectSelection;
0117     StringObjectFunction<std::vector<TOutputCandidateType> > m_combinedObjectSortFunction;
0118     std::map<std::string, std::shared_ptr<StringObjectFunction<std::vector<TOutputCandidateType> > > >
0119         m_plottersCombinedObject;
0120     std::map<std::string, std::shared_ptr<StringObjectFunction<TInputCandidateType> > > m_plottersSingleObject;
0121     /// xxx
0122     static const int SingleObjectPlotter = 0;
0123     static const int CombinedObjectPlotter = 1;
0124     std::map<std::string, int> m_plotterType;
0125     std::vector<edm::ParameterSet> m_combinedObjectDrawables;
0126     std::vector<edm::ParameterSet> m_singleObjectDrawables;  // for all single objects passing preselection
0127     bool m_isSetup;
0128     edm::InputTag m_input;
0129     std::map<std::string, edm::EDGetToken> m_tokens;
0130 
0131   public:
0132     HandlerTemplate(const edm::ParameterSet& iConfig, triggerExpression::Data& eventCache)
0133         : BaseHandler(iConfig, eventCache),
0134           m_singleObjectSelection(iConfig.getParameter<std::string>("singleObjectsPreselection")),
0135           m_combinedObjectSelection(iConfig.getParameter<std::string>("combinedObjectSelection")),
0136           m_combinedObjectSortFunction(iConfig.getParameter<std::string>("combinedObjectSortCriteria")) {
0137       std::string type = iConfig.getParameter<std::string>("handlerType");
0138       if (type != "FromHLT") {
0139         m_input = iConfig.getParameter<edm::InputTag>("inputCol");
0140       }
0141 
0142       m_dqmhistolabel = iConfig.getParameter<std::string>("dqmhistolabel");
0143       m_filterPartialName =
0144           iConfig.getParameter<std::string>("partialFilterName");  // std::string find is used to match filter
0145       m_pathPartialName = iConfig.getParameter<std::string>("partialPathName");
0146       m_combinedObjectDimension = iConfig.getParameter<int>("combinedObjectDimension");
0147       m_combinedObjectDrawables = iConfig.getParameter<std::vector<edm::ParameterSet> >("combinedObjectDrawables");
0148       m_singleObjectDrawables = iConfig.getParameter<std::vector<edm::ParameterSet> >("singleObjectDrawables");
0149       m_isSetup = false;
0150     }
0151     ~HandlerTemplate() override = default;
0152     void book(DQMStore::IBooker& booker) override {
0153       if (!m_isSetup) {
0154         booker.setCurrentFolder(m_dirname);
0155         m_isSetup = true;
0156         std::vector<std::vector<edm::ParameterSet>*> todo(2, (std::vector<edm::ParameterSet>*)nullptr);
0157         todo[CombinedObjectPlotter] = &m_combinedObjectDrawables;
0158         todo[SingleObjectPlotter] = &m_singleObjectDrawables;
0159         for (size_t ti = 0; ti < todo.size(); ++ti) {
0160           for (size_t i = 0; i < todo[ti]->size(); ++i) {
0161             std::string histoName = m_dqmhistolabel + "_" + todo[ti]->at(i).template getParameter<std::string>("name");
0162             std::string expression = todo[ti]->at(i).template getParameter<std::string>("expression");
0163             int bins = todo[ti]->at(i).template getParameter<int>("bins");
0164             double rangeLow = todo[ti]->at(i).template getParameter<double>("min");
0165             double rangeHigh = todo[ti]->at(i).template getParameter<double>("max");
0166             m_histos[histoName] = booker.book1D(histoName, histoName, bins, rangeLow, rangeHigh);
0167             m_plotterType[histoName] = ti;
0168             if (ti == CombinedObjectPlotter) {
0169               auto* func = new StringObjectFunction<std::vector<TOutputCandidateType> >(expression);
0170               m_plottersCombinedObject[histoName] =
0171                   std::shared_ptr<StringObjectFunction<std::vector<TOutputCandidateType> > >(func);
0172             } else {
0173               auto* func = new StringObjectFunction<TInputCandidateType>(expression);
0174               m_plottersSingleObject[histoName] = std::shared_ptr<StringObjectFunction<TInputCandidateType> >(func);
0175             }
0176           }
0177         }
0178       }
0179     }
0180     void getAndStoreTokens(edm::ConsumesCollector&& iC) override {
0181       edm::EDGetTokenT<std::vector<TInputCandidateType> > tok = iC.consumes<std::vector<TInputCandidateType> >(m_input);
0182       m_tokens[m_input.encode()] = edm::EDGetToken(tok);
0183     }
0184 
0185     //#############################################################################
0186     // Count objects. To avoid code duplication we do it in a separate template -
0187     //  - partial specialization not easy...:
0188     // http://stackoverflow.com/questions/21182729/specializing-single-method-in-a-big-template-class
0189     //#############################################################################
0190     template <class T>
0191     int count(const edm::Event& iEvent, InputTag& input, StringCutObjectSelector<T>& sel, float weight) {
0192       int ret = 0;
0193       Handle<std::vector<T> > hIn;
0194       iEvent.getByToken(m_tokens[input.encode()], hIn);
0195       if (!hIn.isValid()) {
0196         edm::LogError("FSQDiJetAve") << "product not found: " << input.encode();
0197         return -1;  // return nonsense value
0198       }
0199       for (size_t i = 0; i < hIn->size(); ++i) {
0200         bool preselection = sel(hIn->at(i));
0201         if (preselection) {
0202           fillSingleObjectPlots(hIn->at(i), weight);
0203           ret += 1;
0204         }
0205       }
0206       return ret;
0207     }
0208 
0209     // FIXME (?): code duplication
0210     void fillSingleObjectPlots(const TInputCandidateType& cand, float weight) {
0211       std::map<std::string, MonitorElement*>::iterator it, itE;
0212       it = m_histos.begin();
0213       itE = m_histos.end();
0214       for (; it != itE; ++it) {
0215         if (m_plotterType[it->first] != SingleObjectPlotter)
0216           continue;
0217         float val = (*m_plottersSingleObject[it->first])(cand);
0218         it->second->Fill(val, weight);
0219       }
0220     }
0221     // Notes:
0222     //  - this function (and specialized versions) are responsible for calling
0223     //     fillSingleObjectPlots for all single objects passing the single
0224     //     object preselection criteria
0225     //  - FIXME this function should take only event/ event setup (?)
0226     //  - FIXME responsibility to apply preselection should be elsewhere
0227     //          hard to fix, since we dont want to copy all objects due to
0228     //          performance reasons
0229     //  - note: implementation below working when in/out types are equal
0230     //          in other cases you must provide specialized version (see below)
0231     void getFilteredCands(TInputCandidateType*,
0232                           std::vector<TOutputCandidateType>& cands,
0233                           const edm::Event& iEvent,
0234                           const edm::EventSetup& iSetup,
0235                           const HLTConfigProvider& hltConfig,
0236                           const trigger::TriggerEvent& trgEvent,
0237                           float weight) {
0238       Handle<std::vector<TInputCandidateType> > hIn;
0239       iEvent.getByToken(m_tokens[m_input.encode()], hIn);
0240 
0241       if (!hIn.isValid()) {
0242         edm::LogError("FSQDiJetAve") << "product not found: " << m_input.encode();
0243         return;
0244       }
0245 
0246       for (size_t i = 0; i < hIn->size(); ++i) {
0247         bool preselection = m_singleObjectSelection(hIn->at(i));
0248         if (preselection) {
0249           fillSingleObjectPlots(hIn->at(i), weight);
0250           cands.push_back(hIn->at(i));
0251         }
0252       }
0253     }
0254 
0255     std::vector<std::string> findPathAndFilter(const HLTConfigProvider& hltConfig) {
0256       std::vector<std::string> ret(2, "");
0257       std::string filterFullName = "";
0258       std::string pathFullName = "";
0259       std::vector<std::string> filtersForThisPath;
0260       //int pathIndex = -1;
0261       int numPathMatches = 0;
0262       int numFilterMatches = 0;
0263       for (size_t i = 0; i < hltConfig.size(); ++i) {
0264         if (hltConfig.triggerName(i).find(m_pathPartialName) == std::string::npos)
0265           continue;
0266         pathFullName = hltConfig.triggerName(i);
0267         //pathIndex = i;
0268         ++numPathMatches;
0269         std::vector<std::string> moduleLabels = hltConfig.moduleLabels(i);
0270         for (auto& moduleLabel : moduleLabels) {
0271           if ("EDFilter" == hltConfig.moduleEDMType(moduleLabel)) {
0272             filtersForThisPath.push_back(moduleLabel);
0273             if (moduleLabel.find(m_filterPartialName) != std::string::npos) {
0274               filterFullName = moduleLabel;
0275               ++numFilterMatches;
0276             }
0277           }
0278         }
0279       }
0280 
0281       // LogWarning or LogError?
0282       if (numPathMatches != 1) {
0283         edm::LogInfo("FSQDiJetAve") << "Problem: found " << numPathMatches << " paths matching " << m_pathPartialName
0284                                     << std::endl;
0285         return ret;
0286       }
0287       ret[0] = pathFullName;
0288       if (numFilterMatches != 1) {
0289         edm::LogError("FSQDiJetAve") << "Problem: found " << numFilterMatches << " filter matching "
0290                                      << m_filterPartialName << " in path " << m_pathPartialName << std::endl;
0291         return ret;
0292       }
0293       ret[1] = filterFullName;
0294       return ret;
0295     }
0296 
0297     void analyze(const edm::Event& iEvent,
0298                  const edm::EventSetup& iSetup,
0299                  const HLTConfigProvider& hltConfig,
0300                  const trigger::TriggerEvent& trgEvent,
0301                  const edm::TriggerResults& triggerResults,
0302                  const edm::TriggerNames& triggerNames,
0303                  float weight) override {
0304       size_t found = 0;
0305       for (size_t i = 0; i < triggerNames.size(); ++i) {
0306         auto itUsedPaths = m_usedPaths.begin();
0307         for (; itUsedPaths != m_usedPaths.end(); ++itUsedPaths) {
0308           if (triggerNames.triggerName(i).find(*itUsedPaths) != std::string::npos) {
0309             ++found;
0310             break;
0311           }
0312         }
0313 
0314         if (found == m_usedPaths.size())
0315           break;
0316       }
0317       if (found != m_usedPaths.size()) {
0318         edm::LogInfo("FSQDiJetAve") << "One of requested paths not found, skipping event";
0319         return;
0320       }
0321       if (m_eventCache->configurationUpdated()) {
0322         m_expression->init(*m_eventCache);
0323       }
0324       if (not(*m_expression)(*m_eventCache))
0325         return;
0326 
0327       /*
0328             std::vector<std::string> pathAndFilter = findPathAndFilter(hltConfig);
0329 
0330             std::string pathFullName = pathAndFilter[0];
0331             if (pathFullName == "") {
0332                 return;
0333             }
0334             unsigned indexNum = triggerNames.triggerIndex(pathFullName);
0335             if(indexNum >= triggerNames.size()){
0336                   edm::LogError("FSQDiJetAve") << "Problem determining trigger index for " << pathFullName << " " << m_pathPartialName;
0337             }
0338             if (!triggerResults.accept(indexNum)) return;*/
0339 
0340       std::vector<TOutputCandidateType> cands;
0341       getFilteredCands((TInputCandidateType*)nullptr, cands, iEvent, iSetup, hltConfig, trgEvent, weight);
0342 
0343       if (cands.empty())
0344         return;
0345 
0346       std::vector<TOutputCandidateType> bestCombinationFromCands = getBestCombination(cands);
0347       if (bestCombinationFromCands.empty())
0348         return;
0349 
0350       // plot
0351       std::map<std::string, MonitorElement*>::iterator it, itE;
0352       it = m_histos.begin();
0353       itE = m_histos.end();
0354       for (; it != itE; ++it) {
0355         if (m_plotterType[it->first] != CombinedObjectPlotter)
0356           continue;
0357         float val = (*m_plottersCombinedObject[it->first])(bestCombinationFromCands);
0358         it->second->Fill(val, weight);
0359       }
0360     }
0361 
0362     std::vector<TOutputCandidateType> getBestCombination(std::vector<TOutputCandidateType>& cands) {
0363       int columnSize = cands.size();
0364       std::vector<int> currentCombination(m_combinedObjectDimension, 0);
0365       std::vector<int> bestCombination(m_combinedObjectDimension, -1);
0366 
0367       int maxCombinations = 1;
0368       int cnt = 0;
0369       while (cnt < m_combinedObjectDimension) {
0370         cnt += 1;
0371         maxCombinations *= columnSize;
0372       }
0373 
0374       cnt = 0;
0375       float bestCombinedCandVal = -1;
0376       while (cnt < maxCombinations) {
0377         cnt += 1;
0378 
0379         // 1. Check if current combination contains duplicates
0380         std::vector<int> currentCombinationCopy(currentCombination);
0381         std::vector<int>::iterator it;
0382         std::sort(currentCombinationCopy.begin(), currentCombinationCopy.end());
0383         it = std::unique(currentCombinationCopy.begin(), currentCombinationCopy.end());
0384         currentCombinationCopy.resize(std::distance(currentCombinationCopy.begin(), it));
0385         bool duplicatesPresent = currentCombination.size() != currentCombinationCopy.size();
0386         // 2. If no duplicates found -
0387         //          - check if current combination passes the cut
0388         //          - rank current combination
0389         if (!duplicatesPresent) {  // no duplicates, we can consider this combined object
0390           std::vector<TOutputCandidateType> currentCombinationFromCands;
0391           currentCombinationFromCands.reserve(m_combinedObjectDimension);
0392           for (int i = 0; i < m_combinedObjectDimension; ++i) {
0393             currentCombinationFromCands.push_back(cands.at(currentCombination.at(i)));
0394           }
0395           bool isOK = m_combinedObjectSelection(currentCombinationFromCands);
0396           if (isOK) {
0397             float curVal = m_combinedObjectSortFunction(currentCombinationFromCands);
0398             // FIXME
0399             if (curVal < 0) {
0400               edm::LogError("FSQDiJetAve")
0401                   << "Problem: ranking function returned negative value: " << curVal << std::endl;
0402             } else if (curVal > bestCombinedCandVal) {
0403               //std::cout << curVal << " " << bestCombinedCandVal << std::endl;
0404               bestCombinedCandVal = curVal;
0405               bestCombination = currentCombination;
0406             }
0407           }
0408         }
0409         // 3. Prepare next combination to test
0410         //    note to future self: less error prone method with modulo
0411         currentCombination.at(m_combinedObjectDimension - 1) += 1;  // increase last number
0412         int carry = 0;
0413         for (int i = m_combinedObjectDimension - 1; i >= 0;
0414              --i) {  // iterate over all numbers, check if we are out of range
0415           currentCombination.at(i) += carry;
0416           carry = 0;
0417           if (currentCombination.at(i) >= columnSize) {
0418             carry = 1;
0419             currentCombination.at(i) = 0;
0420           }
0421         }
0422       }  // combinations loop ends
0423 
0424       std::vector<TOutputCandidateType> bestCombinationFromCands;
0425       if (!bestCombination.empty() && bestCombination.at(0) >= 0) {
0426         for (int i = 0; i < m_combinedObjectDimension; ++i) {
0427           bestCombinationFromCands.push_back(cands.at(bestCombination.at(i)));
0428         }
0429       }
0430       return bestCombinationFromCands;
0431     }
0432   };
0433   //#############################################################################
0434   // Read any object inheriting from reco::Candidate. Save p4
0435   //
0436   //  problem: for reco::Candidate there is no reflex dictionary, so selector
0437   //  wont work
0438   //#############################################################################
0439   template <>
0440   void HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector>::getAndStoreTokens(
0441       edm::ConsumesCollector&& iC) {
0442     edm::EDGetTokenT<View<reco::Candidate> > tok = iC.consumes<View<reco::Candidate> >(m_input);
0443     m_tokens[m_input.encode()] = edm::EDGetToken(tok);
0444   }
0445   template <>
0446   void HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector>::getFilteredCands(
0447       reco::Candidate::LorentzVector*,  // pass a dummy pointer, makes possible to select correct getFilteredCands
0448       std::vector<reco::Candidate::LorentzVector>& cands,  // output collection
0449       const edm::Event& iEvent,
0450       const edm::EventSetup& iSetup,
0451       const HLTConfigProvider& hltConfig,
0452       const trigger::TriggerEvent& trgEvent,
0453       float weight) {
0454     Handle<View<reco::Candidate> > hIn;
0455     iEvent.getByToken(m_tokens[m_input.encode()], hIn);
0456     if (!hIn.isValid()) {
0457       edm::LogError("FSQDiJetAve") << "product not found: " << m_input.encode();
0458       return;
0459     }
0460     for (auto const& i : *hIn) {
0461       bool preselection = m_singleObjectSelection(i.p4());
0462       if (preselection) {
0463         fillSingleObjectPlots(i.p4(), weight);
0464         cands.push_back(i.p4());
0465       }
0466     }
0467   }
0468   //#############################################################################
0469   //
0470   // Count any object inheriting from reco::Track. Save into std::vector<int>
0471   // note: this is similar to recoCand counter (code duplication is hard to
0472   //       avoid in this case)
0473   //
0474   //#############################################################################
0475   template <>
0476   void HandlerTemplate<reco::Track, int>::getFilteredCands(
0477       reco::Track*,             // pass a dummy pointer, makes possible to select correct getFilteredCands
0478       std::vector<int>& cands,  // output collection
0479       const edm::Event& iEvent,
0480       const edm::EventSetup& iSetup,
0481       const HLTConfigProvider& hltConfig,
0482       const trigger::TriggerEvent& trgEvent,
0483       float weight) {
0484     cands.clear();
0485     cands.push_back(count<reco::Track>(iEvent, m_input, m_singleObjectSelection, weight));
0486   }
0487   template <>
0488   void HandlerTemplate<reco::GenParticle, int>::getFilteredCands(
0489       reco::GenParticle*,       // pass a dummy pointer, makes possible to select correct getFilteredCands
0490       std::vector<int>& cands,  // output collection
0491       const edm::Event& iEvent,
0492       const edm::EventSetup& iSetup,
0493       const HLTConfigProvider& hltConfig,
0494       const trigger::TriggerEvent& trgEvent,
0495       float weight) {
0496     cands.clear();
0497     cands.push_back(count<reco::GenParticle>(iEvent, m_input, m_singleObjectSelection, weight));
0498   }
0499   //#############################################################################
0500   //
0501   // Count any object inheriting from reco::Track that is not to distant from
0502   // selected vertex. Save into std::vector<int>
0503   // note: this is similar to recoCand counter (code duplication is hard to
0504   //       avoid in this case)
0505   //
0506   //#############################################################################
0507   template <>
0508   void HandlerTemplate<reco::Track, int, BestVertexMatching>::getAndStoreTokens(edm::ConsumesCollector&& iC) {
0509     edm::EDGetTokenT<std::vector<reco::Track> > tok = iC.consumes<std::vector<reco::Track> >(m_input);
0510     m_tokens[m_input.encode()] = edm::EDGetToken(tok);
0511 
0512     edm::InputTag lVerticesTag = m_pset.getParameter<edm::InputTag>("vtxCollection");
0513     edm::EDGetTokenT<reco::VertexCollection> tok2 = iC.consumes<reco::VertexCollection>(lVerticesTag);
0514     m_tokens[lVerticesTag.encode()] = edm::EDGetToken(tok2);
0515   }
0516 
0517   template <>
0518   void HandlerTemplate<reco::Track, int, BestVertexMatching>::getFilteredCands(
0519       reco::Track*,             // pass a dummy pointer, makes possible to select correct getFilteredCands
0520       std::vector<int>& cands,  // output collection
0521       const edm::Event& iEvent,
0522       const edm::EventSetup& iSetup,
0523       const HLTConfigProvider& hltConfig,
0524       const trigger::TriggerEvent& trgEvent,
0525       float weight) {
0526     // this is not elegant, but should be thread safe
0527     static const edm::InputTag lVerticesTag = m_pset.getParameter<edm::InputTag>("vtxCollection");
0528     static const int lMinNDOF = m_pset.getParameter<int>("minNDOF");                        //7
0529     static const double lMaxZ = m_pset.getParameter<double>("maxZ");                        // 15
0530     static const double lMaxDZ = m_pset.getParameter<double>("maxDZ");                      // 0.12
0531     static const double lMaxDZ2dzsigma = m_pset.getParameter<double>("maxDZ2dzsigma");      // 3
0532     static const double lMaxDXY = m_pset.getParameter<double>("maxDXY");                    // 0.12
0533     static const double lMaxDXY2dxysigma = m_pset.getParameter<double>("maxDXY2dxysigma");  // 3
0534 
0535     cands.clear();
0536     cands.push_back(0);
0537 
0538     edm::Handle<reco::VertexCollection> vertices;
0539     iEvent.getByToken(m_tokens[lVerticesTag.encode()], vertices);
0540 
0541     //double bestvz=-999.9, bestvx=-999.9, bestvy=-999.9;
0542 
0543     double dxy, dz, dzsigma, dxysigma;
0544     math::XYZPoint vtxPoint(0.0, 0.0, 0.0);
0545     double vzErr = 0.0, vxErr = 0.0, vyErr = 0.0;
0546 
0547     // take first vertex passing the criteria
0548     int bestVtx = -1;
0549     for (size_t i = 0; i < vertices->size(); ++i) {
0550       if (vertices->at(i).ndof() < lMinNDOF)
0551         continue;
0552       if (fabs(vertices->at(i).z()) > lMaxZ)
0553         continue;
0554 
0555       vtxPoint = vertices->at(i).position();
0556       vzErr = vertices->at(i).zError();
0557       vxErr = vertices->at(i).xError();
0558       vyErr = vertices->at(i).yError();
0559       bestVtx = i;
0560       break;
0561     }
0562     if (bestVtx < 0)
0563       return;
0564     // const reco::Vertex & vtx = vertices->at(bestVtx);
0565 
0566     Handle<std::vector<reco::Track> > hIn;
0567     iEvent.getByToken(m_tokens[m_input.encode()], hIn);
0568     if (!hIn.isValid()) {
0569       edm::LogError("FSQDiJetAve") << "product not found: " << m_input.encode();
0570       return;
0571     }
0572 
0573     for (auto const& i : *hIn) {
0574       if (!m_singleObjectSelection(i))
0575         continue;
0576       dxy = 0.0, dz = 0.0, dxysigma = 0.0, dzsigma = 0.0;
0577       dxy = -1. * i.dxy(vtxPoint);
0578       dz = i.dz(vtxPoint);
0579       dxysigma = sqrt(i.dxyError() * i.dxyError() + vxErr * vyErr);
0580       dzsigma = sqrt(i.dzError() * i.dzError() + vzErr * vzErr);
0581 
0582       if (fabs(dz) > lMaxDZ)
0583         continue;  // TODO...
0584       if (fabs(dz / dzsigma) > lMaxDZ2dzsigma)
0585         continue;
0586       if (fabs(dxy) > lMaxDXY)
0587         continue;
0588       if (fabs(dxy / dxysigma) > lMaxDXY2dxysigma)
0589         continue;
0590 
0591       cands.at(0) += 1;
0592     }  //loop over tracks
0593   }
0594   //#############################################################################
0595   //
0596   // Apply JEC to PFJets
0597   //
0598   //#############################################################################
0599   template <>
0600   void HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC>::getAndStoreTokens(edm::ConsumesCollector&& iC) {
0601     edm::EDGetTokenT<std::vector<reco::PFJet> > tok = iC.consumes<std::vector<reco::PFJet> >(m_input);
0602     m_tokens[m_input.encode()] = edm::EDGetToken(tok);
0603 
0604     edm::InputTag jetCorTag = m_pset.getParameter<edm::InputTag>("PFJetCorLabel");
0605     edm::EDGetTokenT<reco::JetCorrector> jetcortoken = iC.consumes<reco::JetCorrector>(jetCorTag);
0606     m_tokens[jetCorTag.encode()] = edm::EDGetToken(jetcortoken);
0607   }
0608 
0609   template <>
0610   void HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC>::getFilteredCands(
0611       reco::PFJet*,                     // pass a dummy pointer, makes possible to select correct getFilteredCands
0612       std::vector<reco::PFJet>& cands,  // output collection
0613       const edm::Event& iEvent,
0614       const edm::EventSetup& iSetup,
0615       const HLTConfigProvider& hltConfig,
0616       const trigger::TriggerEvent& trgEvent,
0617       float weight) {
0618     cands.clear();
0619     static const edm::InputTag jetCorTag = m_pset.getParameter<edm::InputTag>("PFJetCorLabel");
0620     edm::Handle<reco::JetCorrector> pfcorrector;
0621     iEvent.getByToken(m_tokens[jetCorTag.encode()], pfcorrector);
0622 
0623     Handle<std::vector<reco::PFJet> > hIn;
0624     iEvent.getByToken(m_tokens[m_input.encode()], hIn);
0625 
0626     if (!hIn.isValid()) {
0627       edm::LogError("FSQDiJetAve") << "product not found: " << m_input.encode();
0628       return;
0629     }
0630 
0631     for (auto const& i : *hIn) {
0632       double scale = pfcorrector->correction(i);
0633       reco::PFJet newPFJet(scale * i.p4(), i.vertex(), i.getSpecific(), i.getJetConstituents());
0634 
0635       bool preselection = m_singleObjectSelection(newPFJet);
0636       if (preselection) {
0637         fillSingleObjectPlots(newPFJet, weight);
0638         cands.push_back(newPFJet);
0639       }
0640     }
0641   }
0642   //#############################################################################
0643   //
0644   // Count any object inheriting from reco::Candidate. Save into std::vector<int>
0645   //  same problem as for reco::Candidate handler ()
0646   //
0647   //#############################################################################
0648   template <>
0649   void HandlerTemplate<reco::Candidate::LorentzVector, int>::getAndStoreTokens(edm::ConsumesCollector&& iC) {
0650     edm::EDGetTokenT<View<reco::Candidate> > tok = iC.consumes<View<reco::Candidate> >(m_input);
0651     m_tokens[m_input.encode()] = edm::EDGetToken(tok);
0652   }
0653   template <>
0654   void HandlerTemplate<reco::Candidate::LorentzVector, int>::getFilteredCands(
0655       reco::Candidate::LorentzVector*,  // pass a dummy pointer, makes possible to select correct getFilteredCands
0656       std::vector<int>& cands,          // output collection
0657       const edm::Event& iEvent,
0658       const edm::EventSetup& iSetup,
0659       const HLTConfigProvider& hltConfig,
0660       const trigger::TriggerEvent& trgEvent,
0661       float weight) {
0662     cands.clear();
0663     cands.push_back(0);
0664 
0665     Handle<View<reco::Candidate> > hIn;
0666     iEvent.getByToken(m_tokens[m_input.encode()], hIn);
0667     if (!hIn.isValid()) {
0668       edm::LogError("FSQDiJetAve") << "product not found: " << m_input.encode();
0669       return;
0670     }
0671     for (auto const& i : *hIn) {
0672       bool preselection = m_singleObjectSelection(i.p4());
0673       if (preselection) {
0674         fillSingleObjectPlots(i.p4(), weight);
0675         cands.at(0) += 1;
0676       }
0677     }
0678   }
0679   //#############################################################################
0680   //
0681   // Read and save trigger::TriggerObject from triggerEvent
0682   //
0683   //#############################################################################
0684   template <>
0685   void HandlerTemplate<trigger::TriggerObject, trigger::TriggerObject>::getFilteredCands(
0686       trigger::TriggerObject*,
0687       std::vector<trigger::TriggerObject>& cands,
0688       const edm::Event& iEvent,
0689       const edm::EventSetup& iSetup,
0690       const HLTConfigProvider& hltConfig,
0691       const trigger::TriggerEvent& trgEvent,
0692       float weight) {
0693     // 1. Find matching path. Inside matchin path find matching filter
0694     std::string filterFullName = findPathAndFilter(hltConfig)[1];
0695     if (filterFullName.empty()) {
0696       return;
0697     }
0698 
0699     // 2. Fetch HLT objects saved by selected filter. Save those fullfilling preselection
0700     //      objects are saved in cands variable
0701     const std::string& process = trgEvent.usedProcessName();  // broken?
0702     edm::InputTag hltTag(filterFullName, "", process);
0703 
0704     const int hltIndex = trgEvent.filterIndex(hltTag);
0705     if (hltIndex >= trgEvent.sizeFilters()) {
0706       edm::LogInfo("FSQDiJetAve") << "Cannot determine hlt index for |" << filterFullName << "|" << process;
0707       return;
0708     }
0709 
0710     const trigger::TriggerObjectCollection& toc(trgEvent.getObjects());
0711     const trigger::Keys& khlt = trgEvent.filterKeys(hltIndex);
0712 
0713     auto kj = khlt.begin();
0714 
0715     for (; kj != khlt.end(); ++kj) {
0716       bool preselection = m_singleObjectSelection(toc[*kj]);
0717       if (preselection) {
0718         fillSingleObjectPlots(toc[*kj], weight);
0719         cands.push_back(toc[*kj]);
0720       }
0721     }
0722   }
0723 
0724   typedef HandlerTemplate<trigger::TriggerObject, trigger::TriggerObject> HLTHandler;
0725   typedef HandlerTemplate<reco::Candidate::LorentzVector, reco::Candidate::LorentzVector>
0726       RecoCandidateHandler;  // in fact reco::Candidate, reco::Candidate::LorentzVector
0727   typedef HandlerTemplate<reco::PFJet, reco::PFJet> RecoPFJetHandler;
0728   typedef HandlerTemplate<reco::PFJet, reco::PFJet, ApplyJEC> RecoPFJetWithJECHandler;
0729   typedef HandlerTemplate<reco::Track, reco::Track> RecoTrackHandler;
0730   typedef HandlerTemplate<reco::Photon, reco::Photon> RecoPhotonHandler;
0731   typedef HandlerTemplate<reco::Muon, reco::Muon> RecoMuonHandler;
0732   typedef HandlerTemplate<reco::GenParticle, reco::GenParticle> RecoGenParticleHandler;
0733   typedef HandlerTemplate<reco::Candidate::LorentzVector, int> RecoCandidateCounter;
0734   typedef HandlerTemplate<reco::Track, int> RecoTrackCounter;
0735   typedef HandlerTemplate<reco::Track, int, BestVertexMatching> RecoTrackCounterWithVertexConstraint;
0736   typedef HandlerTemplate<reco::GenParticle, int> RecoGenParticleCounter;
0737 }  // namespace FSQ
0738 //################################################################################################
0739 //
0740 // Plugin functions
0741 //
0742 //################################################################################################
0743 FSQDiJetAve::FSQDiJetAve(const edm::ParameterSet& iConfig)
0744     : m_eventCache(iConfig.getParameterSet("triggerConfiguration"), consumesCollector()) {
0745   m_useGenWeight = iConfig.getParameter<bool>("useGenWeight");
0746   if (m_useGenWeight) {
0747     m_genEvInfoToken = consumes<GenEventInfoProduct>(edm::InputTag("generator"));
0748   }
0749 
0750   triggerSummaryLabel_ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
0751   triggerResultsLabel_ = iConfig.getParameter<edm::InputTag>("triggerResultsLabel");
0752   triggerSummaryToken = consumes<trigger::TriggerEvent>(triggerSummaryLabel_);
0753   triggerResultsToken = consumes<edm::TriggerResults>(triggerResultsLabel_);
0754 
0755   triggerSummaryFUToken = consumes<trigger::TriggerEvent>(
0756       edm::InputTag(triggerSummaryLabel_.label(), triggerSummaryLabel_.instance(), std::string("FU")));
0757   triggerResultsFUToken = consumes<edm::TriggerResults>(
0758       edm::InputTag(triggerResultsLabel_.label(), triggerResultsLabel_.instance(), std::string("FU")));
0759 
0760   std::vector<edm::ParameterSet> todo = iConfig.getParameter<std::vector<edm::ParameterSet> >("todo");
0761   for (const auto& pset : todo) {
0762     std::string type = pset.getParameter<std::string>("handlerType");
0763     if (type == "FromHLT") {
0764       m_handlers.push_back(std::make_shared<FSQ::HLTHandler>(pset, m_eventCache));
0765     } else if (type == "RecoCandidateCounter") {
0766       m_handlers.push_back(std::make_shared<FSQ::RecoCandidateCounter>(pset, m_eventCache));
0767     } else if (type == "RecoTrackCounter") {
0768       m_handlers.push_back(std::make_shared<FSQ::RecoTrackCounter>(pset, m_eventCache));
0769     } else if (type == "RecoTrackCounterWithVertexConstraint") {
0770       m_handlers.push_back(std::make_shared<FSQ::RecoTrackCounterWithVertexConstraint>(pset, m_eventCache));
0771     } else if (type == "FromRecoCandidate") {
0772       m_handlers.push_back(std::make_shared<FSQ::RecoCandidateHandler>(pset, m_eventCache));
0773     } else if (type == "RecoPFJet") {
0774       m_handlers.push_back(std::make_shared<FSQ::RecoPFJetHandler>(pset, m_eventCache));
0775     } else if (type == "RecoPFJetWithJEC") {
0776       m_handlers.push_back(std::make_shared<FSQ::RecoPFJetWithJECHandler>(pset, m_eventCache));
0777     } else if (type == "RecoTrack") {
0778       m_handlers.push_back(std::make_shared<FSQ::RecoTrackHandler>(pset, m_eventCache));
0779     } else if (type == "RecoPhoton") {
0780       m_handlers.push_back(std::make_shared<FSQ::RecoPhotonHandler>(pset, m_eventCache));
0781     } else if (type == "RecoMuon") {
0782       m_handlers.push_back(std::make_shared<FSQ::RecoMuonHandler>(pset, m_eventCache));
0783     } else if (type == "RecoGenParticleCounter") {
0784       m_handlers.push_back(std::make_shared<FSQ::RecoGenParticleCounter>(pset, m_eventCache));
0785     } else if (type == "RecoGenParticleHandler") {
0786       m_handlers.push_back(std::make_shared<FSQ::RecoGenParticleHandler>(pset, m_eventCache));
0787     } else {
0788       throw cms::Exception("FSQ DQM handler not know: " + type);
0789     }
0790   }
0791   for (auto& m_handler : m_handlers) {
0792     m_handler->getAndStoreTokens(consumesCollector());
0793   }
0794 }
0795 
0796 FSQDiJetAve::~FSQDiJetAve() = default;
0797 
0798 void FSQDiJetAve::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0799   using namespace edm;
0800   if (not m_eventCache.setEvent(iEvent, iSetup)) {
0801     edm::LogError("FSQDiJetAve") << "Could not setup the filter";
0802   }
0803 
0804   //---------- triggerResults ----------
0805   iEvent.getByToken(triggerResultsToken, m_triggerResults);
0806   if (!m_triggerResults.isValid()) {
0807     iEvent.getByToken(triggerResultsFUToken, m_triggerResults);
0808     if (!m_triggerResults.isValid()) {
0809       edm::LogError("FSQDiJetAve") << "TriggerResults not valid, skippng event";
0810       return;
0811     }
0812   }
0813 
0814   //---------- triggerResults ----------
0815   if (m_triggerResults.isValid()) {
0816     m_triggerNames = iEvent.triggerNames(*m_triggerResults);
0817   } else {
0818     edm::LogError("FSQDiJetAve") << "TriggerResults not found";
0819     return;
0820   }
0821 
0822   //---------- triggerSummary ----------
0823   iEvent.getByToken(triggerSummaryToken, m_trgEvent);
0824   if (!m_trgEvent.isValid()) {
0825     iEvent.getByToken(triggerSummaryFUToken, m_trgEvent);
0826     if (!m_trgEvent.isValid()) {
0827       edm::LogInfo("FSQDiJetAve") << "TriggerEvent not found, ";
0828       return;
0829     }
0830   }
0831 
0832   float weight = 1.;
0833   if (m_useGenWeight) {
0834     edm::Handle<GenEventInfoProduct> hGW;
0835     iEvent.getByToken(m_genEvInfoToken, hGW);
0836     weight = hGW->weight();
0837   }
0838 
0839   for (auto& m_handler : m_handlers) {
0840     m_handler->analyze(
0841         iEvent, iSetup, m_hltConfig, *m_trgEvent.product(), *m_triggerResults.product(), m_triggerNames, weight);
0842   }
0843 }
0844 // ------------ method called when starting to processes a run  ------------
0845 //*
0846 void FSQDiJetAve::dqmBeginRun(edm::Run const& run, edm::EventSetup const& c) {
0847   bool changed(true);
0848   std::string processName = triggerResultsLabel_.process();
0849   if (m_hltConfig.init(run, c, processName, changed)) {
0850     LogDebug("FSQDiJetAve") << "HLTConfigProvider failed to initialize.";
0851   }
0852 }
0853 void FSQDiJetAve::bookHistograms(DQMStore::IBooker& booker, edm::Run const& run, edm::EventSetup const& c) {
0854   for (auto& m_handler : m_handlers) {
0855     m_handler->book(booker);
0856   }
0857 }
0858 //*/
0859 // ------------ method called when ending the processing of a run  ------------
0860 /*
0861 void 
0862 FSQDiJetAve::endRun(edm::Run const&, edm::EventSetup const&)
0863 {
0864 }
0865 */
0866 
0867 // ------------ method called when starting to processes a luminosity block  ------------
0868 /*
0869 void 
0870 FSQDiJetAve::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0871 {
0872 }
0873 */
0874 
0875 // ------------ method called when ending the processing of a luminosity block  ------------
0876 /*
0877 void 
0878 FSQDiJetAve::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0879 {}
0880 // */
0881 
0882 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0883 void FSQDiJetAve::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0884   //The following says we do not know what parameters are allowed so do no validation
0885   // Please change this to state exactly what you do use, even if it is no parameters
0886   edm::ParameterSetDescription desc;
0887   desc.setUnknown();
0888   descriptions.addDefault(desc);
0889 }