Trigger/plugins/FSQDiJetAve.cc

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