File indexing completed on 2023-03-17 11:09:17
0001
0002
0003
0004
0005
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
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(""));
0089 desc.add<std::vector<double> >("energyLowEdges", {0.0});
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});
0104 desc.add<edm::InputTag>("l1EGCand", edm::InputTag("hltL1IsoRecoEcalCandidate"));
0105 descriptions.add("hltEgammaGenericQuadraticFilter", desc);
0106 }
0107
0108 HLTEgammaGenericQuadraticFilter::~HLTEgammaGenericQuadraticFilter() {}
0109
0110
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
0120 edm::Ref<reco::RecoEcalCandidateCollection> ref;
0121
0122
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);
0135
0136
0137 edm::Handle<reco::RecoEcalCandidateIsolationMap> depMap;
0138 iEvent.getByToken(varToken_, depMap);
0139
0140
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
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
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.;
0172
0173
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
0206 bool accept(n >= ncandcut_);
0207
0208 return accept;
0209 }
0210
0211
0212 #include "FWCore/Framework/interface/MakerMacros.h"
0213 DEFINE_FWK_MODULE(HLTEgammaGenericQuadraticFilter);