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