Egamma/plugins/HLTEgammaGenericQuadraticFilter.cc

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