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