EgammaHLTProducers/plugins/EgammaHLTFilteredObjProducer.h

0001 #ifndef RECOEGAMMA_EGAMMAHLTPRODUCERS_EGAMMAHLTFILTEREDOBJPRODUCER_H
0002 #define RECOEGAMMA_EGAMMAHLTPRODUCERS_EGAMMAHLTFILTEREDOBJPRODUCER_H
0003
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/stream/EDProducer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013
0014 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
0016 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0017
0018 namespace {
0019   template <typename T>
0020   edm::Handle<T> getHandle(const edm::Event& event, const edm::EDGetTokenT<T>& token) {
0021     edm::Handle<T> handle;
0022     event.getByToken(token, handle);
0023     return handle;
0024   }
0025
0026 }  // namespace
0027
0028 template <typename OutCollType>
0029 class EgammaHLTFilteredObjProducer : public edm::stream::EDProducer<> {
0030 public:
0031   class SelectionCut {
0032   public:
0033     SelectionCut(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0034         : ebCut_(pset.getParameter<edm::ParameterSet>("barrelCut")),
0035           eeCut_(pset.getParameter<edm::ParameterSet>("endcapCut")),
0036           varToken_(iC.consumes<reco::RecoEcalCandidateIsolationMap>(pset.getParameter<edm::InputTag>("var"))) {}
0037
0038     ~SelectionCut() = default;
0039
0040     bool operator()(const reco::RecoEcalCandidateRef& cand) const {
0041       CutValues cut = std::abs(cand->eta()) < kEBEEEtaBound_ ? ebCut_ : eeCut_;
0042       return cut(*cand, getVar(cand));
0043     }
0044
0045     float getVar(const reco::RecoEcalCandidateRef& cand) const {
0046       auto res = varHandle_->find(cand);
0047       if (res != varHandle_->end())
0048         return res->val;
0049       else {
0050         //FIX ME: add some provenance info to this
0051         throw cms::Exception("LogicError") << " candidate not found in collection ";
0052       }
0053     }
0054
0055     void getHandles(const edm::Event& event) { event.getByToken(varToken_, varHandle_); }
0056
0057   private:
0058     struct CutValues {
0059       float cut;
0060       float cutOverE;
0061       float cutOverE2;
0062       bool useEt;
0063       bool doAnd;
0064       std::function<bool(float, float)> compFunc;
0065
0066       CutValues(const edm::ParameterSet& pset)
0067           : cut(pset.getParameter<double>("cut")),
0068             cutOverE(pset.getParameter<double>("cutOverE")),
0069             cutOverE2(pset.getParameter<double>("cutOverE2")),
0070             useEt(pset.getParameter<bool>("useEt")),
0071             doAnd(pset.getParameter<bool>("doAnd")),
0072             compFunc(std::less<float>()) {}
0073
0074       bool operator()(const reco::RecoEcalCandidate& cand, float value) const {
0075         if (!doAnd) {
0076           if (compFunc(value, cut))
0077             return true;
0078           else {
0079             float energyInv = useEt ? 1. / cand.et() : 1. / cand.energy();
0080             return compFunc(value * energyInv, cutOverE) || compFunc(value * energyInv * energyInv, cutOverE2);
0081           }
0082         } else {
0083           float energy = useEt ? cand.et() : cand.energy();
0084           return compFunc(value, cut + cutOverE * energy + cutOverE2 * energy * energy);
0085         }
0086       }
0087     };
0088
0089     CutValues ebCut_;
0090     CutValues eeCut_;
0091     edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> varToken_;
0092     edm::Handle<reco::RecoEcalCandidateIsolationMap> varHandle_;
0093   };
0094
0095   explicit EgammaHLTFilteredObjProducer(const edm::ParameterSet& pset);
0096   ~EgammaHLTFilteredObjProducer() override = default;
0097   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0098
0099   void produce(edm::Event&, const edm::EventSetup&) override;
0100
0101 private:
0102   //different collection types may need to specialise this function
0103   //eg if you wanted superclusters it would output.push_back(cand->superCluster())
0104   static void addObj(const reco::RecoEcalCandidateRef& cand, OutCollType& output) { output.push_back(cand); }
0105
0106   edm::EDGetTokenT<reco::RecoEcalCandidateCollection> candsToken_;
0107   std::vector<SelectionCut> cuts_;
0108   float minEtCutEB_;
0109   float minEtCutEE_;
0110   static constexpr float kEBEEEtaBound_ = 1.479;
0111 };
0112
0113 template <typename OutCollType>
0114 EgammaHLTFilteredObjProducer<OutCollType>::EgammaHLTFilteredObjProducer(const edm::ParameterSet& pset)
0115     : candsToken_(consumes<reco::RecoEcalCandidateCollection>(pset.getParameter<edm::InputTag>("cands"))),
0116       minEtCutEB_(pset.getParameter<double>("minEtCutEB")),
0117       minEtCutEE_(pset.getParameter<double>("minEtCutEE")) {
0118   const auto& cutPsets = pset.getParameter<std::vector<edm::ParameterSet> >("cuts");
0119   for (auto& cutPset : cutPsets) {
0120     cuts_.push_back(SelectionCut(cutPset, consumesCollector()));
0121   }
0122
0123   produces<OutCollType>();
0124 }
0125
0126 template <typename OutCollType>
0127 void EgammaHLTFilteredObjProducer<OutCollType>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0128   edm::ParameterSetDescription desc;
0129   desc.add<edm::InputTag>("cands", edm::InputTag("hltEgammaCandidates"));
0130   desc.add<double>("minEtCutEB", 0);
0131   desc.add<double>("minEtCutEE", 0);
0132
0133   edm::ParameterSetDescription cutsDesc;
0134   edm::ParameterSetDescription regionCutsDesc;
0135   regionCutsDesc.add<double>("cut", -1);
0136   regionCutsDesc.add<double>("cutOverE", -1);
0137   regionCutsDesc.add<double>("cutOverE2", -1);
0138   regionCutsDesc.add<bool>("useEt", false);
0139   regionCutsDesc.add<bool>("doAnd", false);
0140   edm::ParameterSet cutDefaults;
0141   cutDefaults.addParameter<double>("cutOverE", 0.2);
0142   cutDefaults.addParameter<double>("useEt", false);
0143   cutDefaults.addParameter<double>("doAnd", false);
0144
0145   cutsDesc.add<edm::ParameterSetDescription>("barrelCut", regionCutsDesc);
0146   cutsDesc.add<edm::ParameterSetDescription>("endcapCut", regionCutsDesc);
0147   cutsDesc.add<edm::InputTag>("var", edm::InputTag("hltEgammaHoverE"));
0148
0149   edm::ParameterSet defaults;
0150   defaults.addParameter<edm::InputTag>("var", edm::InputTag("hltEgammaHoverE"));
0151   defaults.addParameter<edm::ParameterSet>("barrelCut", cutDefaults);
0152   defaults.addParameter<edm::ParameterSet>("endcapCut", cutDefaults);
0153   desc.addVPSet("cuts", cutsDesc, std::vector<edm::ParameterSet>{defaults});
0154   descriptions.addWithDefaultLabel(desc);
0155 }
0156
0157 template <typename OutCollType>
0158 void EgammaHLTFilteredObjProducer<OutCollType>::produce(edm::Event& iEvent, const edm::EventSetup&) {
0159   for (auto& cut : cuts_)
0160     cut.getHandles(iEvent);
0161   auto candsHandle = getHandle(iEvent, candsToken_);
0162
0163   auto output = std::make_unique<OutCollType>();
0164
0165   for (size_t candNr = 0; candNr < candsHandle->size(); candNr++) {
0166     reco::RecoEcalCandidateRef candRef(candsHandle, candNr);
0167     bool passAllCuts = true;
0168
0169     float minEtCut = std::abs(candRef->eta()) < kEBEEEtaBound_ ? minEtCutEB_ : minEtCutEE_;
0170     if (candRef->et() < minEtCut)
0171       passAllCuts = false;
0172
0173     if (passAllCuts) {
0174       for (const auto& cut : cuts_) {
0175         if (!cut(candRef)) {
0176           passAllCuts = false;
0177           break;
0178         }
0179       }
0180     }
0181     if (passAllCuts)
0182       addObj(candRef, *output);
0183   }
0184
0185   iEvent.put(std::move(output));
0186 }
0187
0188 #endif