File indexing completed on 2024-04-06 12:18:20
0001
0002
0003
0004
0005
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
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(""));
0102 desc.add<std::vector<double> >("energyLowEdges", {0.0});
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});
0126 desc.add<edm::InputTag>("l1EGCand", edm::InputTag("hltL1IsoRecoEcalCandidate"));
0127 descriptions.add("hltEgammaGenericQuadraticEtaFilter", desc);
0128 }
0129
0130 HLTEgammaGenericQuadraticEtaFilter::~HLTEgammaGenericQuadraticEtaFilter() {}
0131
0132
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
0141 edm::Ref<reco::RecoEcalCandidateCollection> ref;
0142
0143
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);
0157
0158
0159 edm::Handle<reco::RecoEcalCandidateIsolationMap> depMap;
0160 iEvent.getByToken(varToken_, depMap);
0161
0162
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
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
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.;
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
0270 bool accept(n >= ncandcut_);
0271
0272 return accept;
0273 }
0274
0275
0276 #include "FWCore/Framework/interface/MakerMacros.h"
0277 DEFINE_FWK_MODULE(HLTEgammaGenericQuadraticEtaFilter);