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