Egamma/plugins/HLTEgammaGenericQuadraticEtaFilter.cc

0001 /** \class HLTEgammaGenericQuadraticEtaFilter
0002  *
0003  *
0004  *  \author Roberto Covarelli (CERN)
0005  *  modified by Chris Tully (Princeton)
0006  */
0007
0008 #include "HLTEgammaGenericQuadraticEtaFilter.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 HLTEgammaGenericQuadraticEtaFilter::HLTEgammaGenericQuadraticEtaFilter(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   etaBoundaryEB12_ = iConfig.getParameter<double>("etaBoundaryEB12");
0036   etaBoundaryEE12_ = iConfig.getParameter<double>("etaBoundaryEE12");
0037
0038   thrRegularEB1_ = iConfig.getParameter<std::vector<double> >("thrRegularEB1");
0039   thrRegularEE1_ = iConfig.getParameter<std::vector<double> >("thrRegularEE1");
0040   thrOverEEB1_ = iConfig.getParameter<std::vector<double> >("thrOverEEB1");
0041   thrOverEEE1_ = iConfig.getParameter<std::vector<double> >("thrOverEEE1");
0042   thrOverE2EB1_ = iConfig.getParameter<std::vector<double> >("thrOverE2EB1");
0043   thrOverE2EE1_ = iConfig.getParameter<std::vector<double> >("thrOverE2EE1");
0044   thrRegularEB2_ = iConfig.getParameter<std::vector<double> >("thrRegularEB2");
0045   thrRegularEE2_ = iConfig.getParameter<std::vector<double> >("thrRegularEE2");
0046   thrOverEEB2_ = iConfig.getParameter<std::vector<double> >("thrOverEEB2");
0047   thrOverEEE2_ = iConfig.getParameter<std::vector<double> >("thrOverEEE2");
0048   thrOverE2EB2_ = iConfig.getParameter<std::vector<double> >("thrOverE2EB2");
0049   thrOverE2EE2_ = iConfig.getParameter<std::vector<double> >("thrOverE2EE2");
0050
0051   ncandcut_ = iConfig.getParameter<int>("ncandcut");
0052
0053   doRhoCorrection_ = iConfig.getParameter<bool>("doRhoCorrection");
0054   rhoMax_ = iConfig.getParameter<double>("rhoMax");
0055   rhoScale_ = iConfig.getParameter<double>("rhoScale");
0056   effectiveAreas_ = iConfig.getParameter<std::vector<double> >("effectiveAreas");
0057   absEtaLowEdges_ = iConfig.getParameter<std::vector<double> >("absEtaLowEdges");
0058
0059   l1EGTag_ = iConfig.getParameter<edm::InputTag>("l1EGCand");
0060
0061   candToken_ = consumes<trigger::TriggerFilterObjectWithRefs>(candTag_);
0062   varToken_ = consumes<reco::RecoEcalCandidateIsolationMap>(varTag_);
0063
0064   if (energyLowEdges_.size() != thrRegularEB1_.size() or energyLowEdges_.size() != thrRegularEE1_.size() or
0065       energyLowEdges_.size() != thrRegularEB2_.size() or energyLowEdges_.size() != thrRegularEE2_.size() or
0066       energyLowEdges_.size() != thrOverEEB1_.size() or energyLowEdges_.size() != thrOverEEE1_.size() or
0067       energyLowEdges_.size() != thrOverEEB2_.size() or energyLowEdges_.size() != thrOverEEE2_.size() or
0068       energyLowEdges_.size() != thrOverE2EB1_.size() or energyLowEdges_.size() != thrOverE2EE1_.size() or
0069       energyLowEdges_.size() != thrOverE2EB2_.size() or energyLowEdges_.size() != thrOverE2EE2_.size())
0070     throw cms::Exception("IncompatibleVects") << "energyLowEdges and threshold vectors should be of the same size. \n";
0071
0072   if (energyLowEdges_.at(0) != 0.0)
0073     throw cms::Exception("IncompleteCoverage") << "energyLowEdges should start from 0. \n";
0074
0075   for (unsigned int aIt = 0; aIt < energyLowEdges_.size() - 1; aIt++) {
0076     if (!(energyLowEdges_.at(aIt) < energyLowEdges_.at(aIt + 1)))
0077       throw cms::Exception("ImproperBinning") << "energyLowEdges entries should be in increasing order. \n";
0078   }
0079
0080   if (doRhoCorrection_) {
0081     rhoToken_ = consumes<double>(rhoTag_);
0082     if (absEtaLowEdges_.size() != effectiveAreas_.size())
0083       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0084
0085     if (absEtaLowEdges_.at(0) != 0.0)
0086       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0087
0088     for (unsigned int bIt = 0; bIt < absEtaLowEdges_.size() - 1; bIt++) {
0089       if (!(absEtaLowEdges_.at(bIt) < absEtaLowEdges_.at(bIt + 1)))
0090         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0091     }
0092   }
0093 }
0094
0095 void HLTEgammaGenericQuadraticEtaFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0096   edm::ParameterSetDescription desc;
0097   makeHLTFilterDescription(desc);
0098   desc.add<edm::InputTag>("candTag", edm::InputTag("hltEGIsolFilter"));
0099   desc.add<edm::InputTag>("varTag", edm::InputTag("hltEGIsol"));
0100   desc.add<edm::InputTag>("rhoTag", edm::InputTag(""));     // No rho correction by default
0101   desc.add<std::vector<double> >("energyLowEdges", {0.0});  // No energy-dependent cuts by default
0102   desc.add<bool>("lessThan", true);
0103   desc.add<bool>("useEt", true);
0104   desc.add<double>("etaBoundaryEB12", 1.0);
0105   desc.add<double>("etaBoundaryEE12", 2.0);
0106   desc.add<std::vector<double> >("thrRegularEB1", {4.0});
0107   desc.add<std::vector<double> >("thrRegularEE1", {6.0});
0108   desc.add<std::vector<double> >("thrOverEEB1", {0.0020});
0109   desc.add<std::vector<double> >("thrOverEEE1", {0.0020});
0110   desc.add<std::vector<double> >("thrOverE2EB1", {0.0});
0111   desc.add<std::vector<double> >("thrOverE2EE1", {0.0});
0112   desc.add<std::vector<double> >("thrRegularEB2", {6.0});
0113   desc.add<std::vector<double> >("thrRegularEE2", {4.0});
0114   desc.add<std::vector<double> >("thrOverEEB2", {0.0020});
0115   desc.add<std::vector<double> >("thrOverEEE2", {0.0020});
0116   desc.add<std::vector<double> >("thrOverE2EB2", {0.0});
0117   desc.add<std::vector<double> >("thrOverE2EE2", {0.0});
0118   desc.add<int>("ncandcut", 1);
0119   desc.add<bool>("doRhoCorrection", false);
0120   desc.add<double>("rhoMax", 9.9999999E7);
0121   desc.add<double>("rhoScale", 1.0);
0122   desc.add<std::vector<double> >("effectiveAreas", {0.0, 0.0});
0123   desc.add<std::vector<double> >("absEtaLowEdges", {0.0, 1.479});  // EB, EE
0124   desc.add<edm::InputTag>("l1EGCand", edm::InputTag("hltL1IsoRecoEcalCandidate"));
0125   descriptions.add("hltEgammaGenericQuadraticEtaFilter", desc);
0126 }
0127
0128 HLTEgammaGenericQuadraticEtaFilter::~HLTEgammaGenericQuadraticEtaFilter() {}
0129
0130 // ------------ method called to produce the data  ------------
0131 bool HLTEgammaGenericQuadraticEtaFilter::hltFilter(edm::Event& iEvent,
0132                                                    const edm::EventSetup& iSetup,
0133                                                    trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0134   using namespace trigger;
0135   if (saveTags()) {
0136     filterproduct.addCollectionTag(l1EGTag_);
0137   }
0138   // Ref to Candidate object to be recorded in filter object
0139   edm::Ref<reco::RecoEcalCandidateCollection> ref;
0140
0141   // Set output format
0142   int trigger_type = trigger::TriggerCluster;
0143   if (saveTags())
0144     trigger_type = trigger::TriggerPhoton;
0145
0146   edm::Handle<trigger::TriggerFilterObjectWithRefs> PrevFilterOutput;
0147
0148   iEvent.getByToken(candToken_, PrevFilterOutput);
0149
0150   std::vector<edm::Ref<reco::RecoEcalCandidateCollection> > recoecalcands;
0151   PrevFilterOutput->getObjects(TriggerCluster, recoecalcands);
0152   if (recoecalcands.empty())
0153     PrevFilterOutput->getObjects(TriggerPhoton,
0154                                  recoecalcands);  //we dont know if its type trigger cluster or trigger photon
0155
0156   //get hold of isolated association map
0157   edm::Handle<reco::RecoEcalCandidateIsolationMap> depMap;
0158   iEvent.getByToken(varToken_, depMap);
0159
0160   // Get rho if needed
0161   edm::Handle<double> rhoHandle;
0162   double rho = 0.0;
0163   if (doRhoCorrection_) {
0164     iEvent.getByToken(rhoToken_, rhoHandle);
0165     rho = *(rhoHandle.product());
0166   }
0167
0168   if (rho > rhoMax_)
0169     rho = rhoMax_;
0170   rho = rho * rhoScale_;
0171
0172   // look at all photons, check cuts and add to filter object
0173   int n = 0;
0174   for (unsigned int i = 0; i < recoecalcands.size(); i++) {
0175     ref = recoecalcands[i];
0176     reco::RecoEcalCandidateIsolationMap::const_iterator mapi = (*depMap).find(ref);
0177
0178     float vali = mapi->val;
0179     float EtaSC = ref->eta();
0180
0181     // Pick the right EA and do rhoCorr
0182     if (doRhoCorrection_) {
0183       auto cIt = std::lower_bound(absEtaLowEdges_.begin(), absEtaLowEdges_.end(), std::abs(EtaSC)) - 1;
0184       vali = vali - (rho * effectiveAreas_.at(std::distance(absEtaLowEdges_.begin(), cIt)));
0185     }
0186
0187     float energy = ref->superCluster()->energy();
0188     if (useEt_)
0189       energy = energy * sin(2 * atan(exp(-EtaSC)));
0190     if (energy < 0.)
0191       energy = 0.; /* first and second order terms assume non-negative energies */
0192
0193     double cutRegularEB1_ = 9999., cutRegularEE1_ = 9999.;
0194     double cutRegularEB2_ = 9999., cutRegularEE2_ = 9999.;
0195     double cutOverEEB1_ = 9999., cutOverEEE1_ = 9999.;
0196     double cutOverEEB2_ = 9999., cutOverEEE2_ = 9999.;
0197     double cutOverE2EB1_ = 9999., cutOverE2EE1_ = 9999.;
0198     double cutOverE2EB2_ = 9999., cutOverE2EE2_ = 9999.;
0199
0200     auto dIt = std::lower_bound(energyLowEdges_.begin(), energyLowEdges_.end(), energy) - 1;
0201     unsigned iEn = std::distance(energyLowEdges_.begin(), dIt);
0202
0203     cutRegularEB1_ = thrRegularEB1_.at(iEn);
0204     cutRegularEB2_ = thrRegularEB2_.at(iEn);
0205     cutRegularEE1_ = thrRegularEE1_.at(iEn);
0206     cutRegularEE2_ = thrRegularEE2_.at(iEn);
0207     cutOverEEB1_ = thrOverEEB1_.at(iEn);
0208     cutOverEEB2_ = thrOverEEB2_.at(iEn);
0209     cutOverEEE1_ = thrOverEEE1_.at(iEn);
0210     cutOverEEE2_ = thrOverEEE2_.at(iEn);
0211     cutOverE2EB1_ = thrOverE2EB1_.at(iEn);
0212     cutOverE2EB2_ = thrOverE2EB2_.at(iEn);
0213     cutOverE2EE1_ = thrOverE2EE1_.at(iEn);
0214     cutOverE2EE2_ = thrOverE2EE2_.at(iEn);
0215
0216     if (lessThan_) {
0217       if (std::abs(EtaSC) < etaBoundaryEB12_) {
0218         if (vali <= cutRegularEB1_ + energy * cutOverEEB1_ + energy * energy * cutOverE2EB1_) {
0219           n++;
0220           filterproduct.addObject(trigger_type, ref);
0221           continue;
0222         }
0223       } else if (std::abs(EtaSC) < 1.479) {
0224         if (vali <= cutRegularEB2_ + energy * cutOverEEB2_ + energy * energy * cutOverE2EB2_) {
0225           n++;
0226           filterproduct.addObject(trigger_type, ref);
0227           continue;
0228         }
0229       } else if (std::abs(EtaSC) < etaBoundaryEE12_) {
0230         if (vali <= cutRegularEE1_ + energy * cutOverEEE1_ + energy * energy * cutOverE2EE1_) {
0231           n++;
0232           filterproduct.addObject(trigger_type, ref);
0233           continue;
0234         }
0235       } else if (vali <= cutRegularEE2_ + energy * cutOverEEE2_ + energy * energy * cutOverE2EE2_) {
0236         n++;
0237         filterproduct.addObject(trigger_type, ref);
0238         continue;
0239       }
0240     } else {
0241       if (std::abs(EtaSC) < etaBoundaryEB12_) {
0242         if (vali >= cutRegularEB1_ + energy * cutOverEEB1_ + energy * energy * cutOverE2EB1_) {
0243           n++;
0244           filterproduct.addObject(trigger_type, ref);
0245           continue;
0246         }
0247       } else if (std::abs(EtaSC) < 1.479) {
0248         if (vali >= cutRegularEB2_ + energy * cutOverEEB2_ + energy * energy * cutOverE2EB2_) {
0249           n++;
0250           filterproduct.addObject(trigger_type, ref);
0251           continue;
0252         }
0253       } else if (std::abs(EtaSC) < etaBoundaryEE12_) {
0254         if (vali >= cutRegularEE1_ + energy * cutOverEEE1_ + energy * energy * cutOverE2EE1_) {
0255           n++;
0256           filterproduct.addObject(trigger_type, ref);
0257           continue;
0258         }
0259       } else if (vali >= cutRegularEE2_ + energy * cutOverEEE2_ + energy * energy * cutOverE2EE2_) {
0260         n++;
0261         filterproduct.addObject(trigger_type, ref);
0262         continue;
0263       }
0264     }
0265   }
0266
0267   // filter decision
0268   bool accept(n >= ncandcut_);
0269
0270   return accept;
0271 }
0272
0273 // declare this class as a framework plugin
0274 #include "FWCore/Framework/interface/MakerMacros.h"
0275 DEFINE_FWK_MODULE(HLTEgammaGenericQuadraticEtaFilter);