Egamma/plugins/HLTGenericFilter.cc

0001 /** \class HLTGenericFilter
0002  *
0003  *
0004  *  \author Roberto Covarelli (CERN)
0005  *
0006  */
0007
0008 #include "HLTGenericFilter.h"
0009
0010 #include "DataFormats/Common/interface/Handle.h"
0011
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0018 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0019 #include "DataFormats/Common/interface/AssociationMap.h"
0020
0021 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0022
0023 //
0024 // constructors and destructor
0025 //
0026 template <typename T1>
0027 HLTGenericFilter<T1>::HLTGenericFilter(const edm::ParameterSet& iConfig) : HLTFilter(iConfig) {
0028   candTag_ = iConfig.template getParameter<edm::InputTag>("candTag");
0029   varTag_ = iConfig.template getParameter<edm::InputTag>("varTag");
0030   l1EGTag_ = iConfig.template getParameter<edm::InputTag>("l1EGCand");
0031   rhoTag_ = iConfig.template getParameter<edm::InputTag>("rhoTag");
0032
0033   energyLowEdges_ = iConfig.template getParameter<std::vector<double>>("energyLowEdges");
0034   lessThan_ = iConfig.template getParameter<bool>("lessThan");
0035   useEt_ = iConfig.template getParameter<bool>("useEt");
0036   useAbs_ = iConfig.template getParameter<bool>("useAbs");
0037   thrRegularEB_ = iConfig.template getParameter<std::vector<double>>("thrRegularEB");
0038   thrRegularEE_ = iConfig.template getParameter<std::vector<double>>("thrRegularEE");
0039   thrOverEEB_ = iConfig.template getParameter<std::vector<double>>("thrOverEEB");
0040   thrOverEEE_ = iConfig.template getParameter<std::vector<double>>("thrOverEEE");
0041   thrOverE2EB_ = iConfig.template getParameter<std::vector<double>>("thrOverE2EB");
0042   thrOverE2EE_ = iConfig.template getParameter<std::vector<double>>("thrOverE2EE");
0043   ncandcut_ = iConfig.template getParameter<int>("ncandcut");
0044
0045   doRhoCorrection_ = iConfig.template getParameter<bool>("doRhoCorrection");
0046   rhoMax_ = iConfig.template getParameter<double>("rhoMax");
0047   rhoScale_ = iConfig.template getParameter<double>("rhoScale");
0048   effectiveAreas_ = iConfig.template getParameter<std::vector<double>>("effectiveAreas");
0049   absEtaLowEdges_ = iConfig.template getParameter<std::vector<double>>("absEtaLowEdges");
0050
0051   candToken_ = consumes<trigger::TriggerFilterObjectWithRefs>(candTag_);
0052   varToken_ = consumes<T1IsolationMap>(varTag_);
0053
0054   if (energyLowEdges_.size() != thrRegularEB_.size() or energyLowEdges_.size() != thrRegularEE_.size() or
0055       energyLowEdges_.size() != thrOverEEB_.size() or energyLowEdges_.size() != thrOverEEE_.size() or
0056       energyLowEdges_.size() != thrOverE2EB_.size() or energyLowEdges_.size() != thrOverE2EE_.size())
0057     throw cms::Exception("IncompatibleVects") << "energyLowEdges and threshold vectors should be of the same size. \n";
0058
0059   if (energyLowEdges_.at(0) != 0.0)
0060     throw cms::Exception("IncompleteCoverage") << "energyLowEdges should start from 0. \n";
0061
0062   for (unsigned int aIt = 0; aIt < energyLowEdges_.size() - 1; aIt++) {
0063     if (!(energyLowEdges_.at(aIt) < energyLowEdges_.at(aIt + 1)))
0064       throw cms::Exception("ImproperBinning") << "energyLowEdges entries should be in increasing order. \n";
0065   }
0066
0067   if (doRhoCorrection_) {
0068     rhoToken_ = consumes<double>(rhoTag_);
0069     if (absEtaLowEdges_.size() != effectiveAreas_.size())
0070       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0071
0072     if (absEtaLowEdges_.at(0) != 0.0)
0073       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0074
0075     for (unsigned int bIt = 0; bIt < absEtaLowEdges_.size() - 1; bIt++) {
0076       if (!(absEtaLowEdges_.at(bIt) < absEtaLowEdges_.at(bIt + 1)))
0077         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0078     }
0079   }
0080 }
0081
0082 template <typename T1>
0083 void HLTGenericFilter<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0084   edm::ParameterSetDescription desc;
0085   makeHLTFilterDescription(desc);
0086   desc.add<edm::InputTag>("candTag", edm::InputTag("hltSingleEgammaEtFilter"));
0087   desc.add<edm::InputTag>("varTag", edm::InputTag("hltSingleEgammaHcalIsol"));
0088   desc.add<edm::InputTag>("rhoTag", edm::InputTag(""));    // No rho correction by default
0089   desc.add<std::vector<double>>("energyLowEdges", {0.0});  // No energy-dependent cuts by default
0090   desc.add<bool>("lessThan", true);
0091   desc.add<bool>("useEt", false);
0092   desc.add<bool>("useAbs", false);
0093   desc.add<std::vector<double>>("thrRegularEB", {0.0});
0094   desc.add<std::vector<double>>("thrRegularEE", {0.0});
0095   desc.add<std::vector<double>>("thrOverEEB", {-1.0});
0096   desc.add<std::vector<double>>("thrOverEEE", {-1.0});
0097   desc.add<std::vector<double>>("thrOverE2EB", {-1.0});
0098   desc.add<std::vector<double>>("thrOverE2EE", {-1.0});
0099   desc.add<int>("ncandcut", 1);
0100   desc.add<bool>("doRhoCorrection", false);
0101   desc.add<double>("rhoMax", 9.9999999E7);
0102   desc.add<double>("rhoScale", 1.0);
0103   desc.add<std::vector<double>>("effectiveAreas", {0.0, 0.0});
0104   desc.add<std::vector<double>>("absEtaLowEdges", {0.0, 1.479});  // EB, EE
0105   desc.add<edm::InputTag>("l1EGCand", edm::InputTag("hltL1IsoRecoEcalCandidate"));
0106   descriptions.add(defaultModuleLabel<HLTGenericFilter<T1>>(), desc);
0107 }
0108
0109 template <typename T1>
0110 HLTGenericFilter<T1>::~HLTGenericFilter() {}
0111
0112 template <typename T1>
0113 float HLTGenericFilter<T1>::getEnergy(T1Ref candRef) const {
0114   return candRef->p();
0115 }
0116
0117 template <>
0118 float HLTGenericFilter<reco::RecoEcalCandidate>::getEnergy(T1Ref candRef) const {
0119   return candRef->superCluster()->energy();
0120 }
0121
0122 template <typename T1>
0123 float HLTGenericFilter<T1>::getEt(T1Ref candRef) const {
0124   return candRef->pt();
0125 }
0126
0127 template <>
0128 float HLTGenericFilter<reco::RecoEcalCandidate>::getEt(T1Ref candRef) const {
0129   return candRef->superCluster()->energy() * sin(2 * atan(exp(-candRef->eta())));
0130 }
0131
0132 // ------------ method called to produce the data  ------------
0133 template <typename T1>
0134 bool HLTGenericFilter<T1>::hltFilter(edm::Event& iEvent,
0135                                      const edm::EventSetup& iSetup,
0136                                      trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0137   using namespace trigger;
0138   if (saveTags()) {
0139     filterproduct.addCollectionTag(l1EGTag_);
0140   }
0141
0142   // Set output format
0143   int trigger_type = trigger::TriggerCluster;
0144   if (saveTags())
0145     trigger_type = trigger::TriggerPhoton;
0146
0147   edm::Handle<trigger::TriggerFilterObjectWithRefs> PrevFilterOutput;
0148   iEvent.getByToken(candToken_, PrevFilterOutput);
0149
0150   std::vector<T1Ref> recoCands;
0151   PrevFilterOutput->getObjects(TriggerCluster, recoCands);
0152   if (recoCands.empty())
0153     PrevFilterOutput->getObjects(TriggerPhoton,
0154                                  recoCands);  //we dont know if its type trigger cluster or trigger photon
0155   if (recoCands.empty()) {
0156     PrevFilterOutput->getObjects(TriggerMuon, recoCands);  //if not a cluster and not a photon then assum it is a muon
0157     trigger_type = trigger::TriggerMuon;
0158   }
0159   //get hold of isolated association map
0160   edm::Handle<T1IsolationMap> depMap;
0161   iEvent.getByToken(varToken_, depMap);
0162
0163   // Get rho if needed
0164   edm::Handle<double> rhoHandle;
0165   double rho = 0.0;
0166   if (doRhoCorrection_) {
0167     iEvent.getByToken(rhoToken_, rhoHandle);
0168     rho = *(rhoHandle.product());
0169   }
0170
0171   if (rho > rhoMax_)
0172     rho = rhoMax_;
0173   rho = rho * rhoScale_;
0174
0175   // look at all photons, check cuts and add to filter object
0176   int n = 0;
0177   for (unsigned int i = 0; i < recoCands.size(); i++) {
0178     // Ref to Candidate object to be recorded in filter object
0179     T1Ref ref = recoCands[i];
0180     typename T1IsolationMap::const_iterator mapi = (*depMap).find(ref);
0181
0182     //should we do the abs before or after the rho corr?
0183     //note: its very unlikely to rho correct a variable that wants to be abs
0184     //decided yes as it seems slightly better this way but can not think of a use case either way
0185     float vali = useAbs_ ? std::abs(mapi->val) : mapi->val;
0186
0187     float EtaSC = ref->eta();
0188
0189     // Pick the right EA and do rhoCorr
0190     if (doRhoCorrection_) {
0191       auto cIt = std::lower_bound(absEtaLowEdges_.begin(), absEtaLowEdges_.end(), std::abs(EtaSC)) - 1;
0192       vali = vali - (rho * effectiveAreas_.at(std::distance(absEtaLowEdges_.begin(), cIt)));
0193     }
0194
0195     float energy;
0196     if (useEt_)
0197       energy = getEt(ref);
0198     else
0199       energy = getEnergy(ref);
0200
0201     // Pick the right cut threshold
0202     double cutRegularEB_ = 9999., cutRegularEE_ = 9999.;
0203     double cutOverEEB_ = 9999., cutOverEEE_ = 9999.;
0204     double cutOverE2EB_ = 9999., cutOverE2EE_ = 9999.;
0205
0206     auto dIt = std::lower_bound(energyLowEdges_.begin(), energyLowEdges_.end(), energy) - 1;
0207     unsigned iEn = std::distance(energyLowEdges_.begin(), dIt);
0208
0209     cutRegularEB_ = thrRegularEB_.at(iEn);
0210     cutRegularEE_ = thrRegularEE_.at(iEn);
0211     cutOverEEB_ = thrOverEEB_.at(iEn);
0212     cutOverEEE_ = thrOverEEE_.at(iEn);
0213     cutOverE2EB_ = thrOverE2EB_.at(iEn);
0214     cutOverE2EE_ = thrOverE2EE_.at(iEn);
0215
0216     if (lessThan_) {
0217       if ((std::abs(EtaSC) < 1.479 && vali <= cutRegularEB_) || (std::abs(EtaSC) >= 1.479 && vali <= cutRegularEE_)) {
0218         n++;
0219         filterproduct.addObject(trigger_type, ref);
0220         continue;
0221       }
0222       if (energy > 0. && (cutOverEEB_ > 0. || cutOverEEE_ > 0. || cutOverE2EB_ > 0. || cutOverE2EE_ > 0.)) {
0223         if ((std::abs(EtaSC) < 1.479 && vali / energy <= cutOverEEB_) ||
0224             (std::abs(EtaSC) >= 1.479 && vali / energy <= cutOverEEE_)) {
0225           n++;
0226           filterproduct.addObject(trigger_type, ref);
0227           continue;
0228         }
0229         if ((std::abs(EtaSC) < 1.479 && vali / (energy * energy) <= cutOverE2EB_) ||
0230             (std::abs(EtaSC) >= 1.479 && vali / (energy * energy) <= cutOverE2EE_)) {
0231           n++;
0232           filterproduct.addObject(trigger_type, ref);
0233         }
0234       }
0235     } else {
0236       if ((std::abs(EtaSC) < 1.479 && vali >= cutRegularEB_) || (std::abs(EtaSC) >= 1.479 && vali >= cutRegularEE_)) {
0237         n++;
0238         filterproduct.addObject(trigger_type, ref);
0239         continue;
0240       }
0241       if (energy > 0. && (cutOverEEB_ > 0. || cutOverEEE_ > 0. || cutOverE2EB_ > 0. || cutOverE2EE_ > 0.)) {
0242         if ((std::abs(EtaSC) < 1.479 && vali / energy >= cutOverEEB_) ||
0243             (std::abs(EtaSC) >= 1.479 && vali / energy >= cutOverEEE_)) {
0244           n++;
0245           filterproduct.addObject(trigger_type, ref);
0246           continue;
0247         }
0248         if ((std::abs(EtaSC) < 1.479 && vali / (energy * energy) >= cutOverE2EB_) ||
0249             (std::abs(EtaSC) >= 1.479 && vali / (energy * energy) >= cutOverE2EE_)) {
0250           n++;
0251           filterproduct.addObject(trigger_type, ref);
0252         }
0253       }
0254     }
0255   }
0256
0257   // filter decision
0258   bool accept(n >= ncandcut_);
0259
0260   return accept;
0261 }
0262
0263 typedef HLTGenericFilter<reco::RecoEcalCandidate> HLTEgammaGenericFilter;
0264 typedef HLTGenericFilter<reco::RecoChargedCandidate> HLTMuonGenericFilter;
0265 DEFINE_FWK_MODULE(HLTEgammaGenericFilter);
0266 DEFINE_FWK_MODULE(HLTMuonGenericFilter);