File indexing completed on 2021-09-02 03:46:02
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0021 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0022 #include "FWCore/Framework/interface/global/EDProducer.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0025 #include "FWCore/Utilities/interface/Exception.h"
0026 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0027 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
0028 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0029 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0030 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0031 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0032 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0033 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0034 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0035
0036 class EgammaHLTHcalVarProducerFromRecHit : public edm::global::EDProducer<> {
0037 public:
0038 explicit EgammaHLTHcalVarProducerFromRecHit(const edm::ParameterSet &);
0039
0040 public:
0041 void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const final;
0042 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0043
0044 private:
0045 const bool doEtSum_;
0046 const EgammaHcalIsolation::arrayHB eThresHB_;
0047 const EgammaHcalIsolation::arrayHB etThresHB_;
0048 const EgammaHcalIsolation::arrayHE eThresHE_;
0049 const EgammaHcalIsolation::arrayHE etThresHE_;
0050 const double innerCone_;
0051 const double outerCone_;
0052 const int depth_;
0053 const int maxSeverityHB_;
0054 const int maxSeverityHE_;
0055 const bool useSingleTower_;
0056 const bool doRhoCorrection_;
0057 const double rhoScale_;
0058 const double rhoMax_;
0059 const std::vector<double> effectiveAreas_;
0060 const std::vector<double> absEtaLowEdges_;
0061 const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateProducer_;
0062 const edm::EDGetTokenT<HBHERecHitCollection> hbheRecHitsTag_;
0063 const edm::EDGetTokenT<double> rhoProducer_;
0064 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0065 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalTopologyToken_;
0066 const edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> hcalChannelQualityToken_;
0067 const edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> hcalSevLvlComputerToken_;
0068 const edm::ESGetToken<CaloTowerConstituentsMap, CaloGeometryRecord> caloTowerConstituentsMapToken_;
0069 const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> putToken_;
0070 };
0071
0072 EgammaHLTHcalVarProducerFromRecHit::EgammaHLTHcalVarProducerFromRecHit(const edm::ParameterSet &config)
0073 : doEtSum_(config.getParameter<bool>("doEtSum")),
0074 eThresHB_(config.getParameter<EgammaHcalIsolation::arrayHB>("eThresHB")),
0075 etThresHB_(config.getParameter<EgammaHcalIsolation::arrayHB>("etThresHB")),
0076 eThresHE_(config.getParameter<EgammaHcalIsolation::arrayHE>("eThresHE")),
0077 etThresHE_(config.getParameter<EgammaHcalIsolation::arrayHE>("etThresHE")),
0078 innerCone_(config.getParameter<double>("innerCone")),
0079 outerCone_(config.getParameter<double>("outerCone")),
0080 depth_(config.getParameter<int>("depth")),
0081 maxSeverityHB_(config.getParameter<int>("maxSeverityHB")),
0082 maxSeverityHE_(config.getParameter<int>("maxSeverityHE")),
0083 useSingleTower_(config.getParameter<bool>("useSingleTower")),
0084 doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0085 rhoScale_(config.getParameter<double>("rhoScale")),
0086 rhoMax_(config.getParameter<double>("rhoMax")),
0087 effectiveAreas_(config.getParameter<std::vector<double> >("effectiveAreas")),
0088 absEtaLowEdges_(config.getParameter<std::vector<double> >("absEtaLowEdges")),
0089 recoEcalCandidateProducer_(consumes(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0090 hbheRecHitsTag_(consumes(config.getParameter<edm::InputTag>("hbheRecHitsTag"))),
0091 rhoProducer_(doRhoCorrection_ ? consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))
0092 : edm::EDGetTokenT<double>()),
0093 caloGeometryToken_{esConsumes()},
0094 hcalTopologyToken_{esConsumes()},
0095 hcalChannelQualityToken_{esConsumes(edm::ESInputTag("", "withTopo"))},
0096 hcalSevLvlComputerToken_{esConsumes()},
0097 caloTowerConstituentsMapToken_{esConsumes()},
0098 putToken_{produces<reco::RecoEcalCandidateIsolationMap>()} {
0099 if (doRhoCorrection_) {
0100 if (absEtaLowEdges_.size() != effectiveAreas_.size()) {
0101 throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0102 }
0103
0104 if (absEtaLowEdges_.at(0) != 0.0) {
0105 throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0106 }
0107
0108 for (unsigned int aIt = 0; aIt < absEtaLowEdges_.size() - 1; aIt++) {
0109 if (!(absEtaLowEdges_.at(aIt) < absEtaLowEdges_.at(aIt + 1))) {
0110 throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0111 }
0112 }
0113 }
0114 }
0115
0116 void EgammaHLTHcalVarProducerFromRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0117 edm::ParameterSetDescription desc;
0118
0119 desc.add<edm::InputTag>("recoEcalCandidateProducer", edm::InputTag("hltRecoEcalCandidate"));
0120 desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0121 desc.add<edm::InputTag>("hbheRecHitsTag", edm::InputTag("hltHbhereco"));
0122 desc.add<bool>("doRhoCorrection", false);
0123 desc.add<double>("rhoMax", 999999.);
0124 desc.add<double>(("rhoScale"), 1.0);
0125
0126 desc.add<std::vector<double> >("eThresHB", {0.1, 0.2, 0.3, 0.3});
0127 desc.add<std::vector<double> >("etThresHB", {0, 0, 0, 0});
0128 desc.add<std::vector<double> >("eThresHE", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0129 desc.add<std::vector<double> >("etThresHE", {0, 0, 0, 0, 0, 0, 0});
0130 desc.add<double>("innerCone", 0);
0131 desc.add<double>("outerCone", 0.14);
0132 desc.add<int>("depth", 0);
0133 desc.add<int>("maxSeverityHB", 9);
0134 desc.add<int>("maxSeverityHE", 9);
0135 desc.add<bool>("doEtSum", false);
0136 desc.add<bool>("useSingleTower", false);
0137 desc.add<std::vector<double> >("effectiveAreas", {0.079, 0.25});
0138 desc.add<std::vector<double> >("absEtaLowEdges", {0.0, 1.479});
0139 descriptions.add("hltEgammaHLTHcalVarProducerFromRecHit", desc);
0140 }
0141
0142 void EgammaHLTHcalVarProducerFromRecHit::produce(edm::StreamID,
0143 edm::Event &iEvent,
0144 const edm::EventSetup &iSetup) const {
0145 auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateProducer_);
0146
0147 double rho = 0.0;
0148
0149 if (doRhoCorrection_) {
0150 rho = iEvent.get(rhoProducer_);
0151 if (rho > rhoMax_) {
0152 rho = rhoMax_;
0153 }
0154 rho = rho * rhoScale_;
0155 }
0156
0157 reco::RecoEcalCandidateIsolationMap isoMap(recoEcalCandHandle);
0158
0159 for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); iRecoEcalCand++) {
0160 reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
0161
0162 float isol = 0;
0163 EgammaHcalIsolation::InclusionRule external;
0164 EgammaHcalIsolation::InclusionRule internal;
0165
0166 if (useSingleTower_) {
0167 if (!doEtSum_) {
0168 external = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0169 internal = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0170 } else {
0171 external = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0172 internal = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0173 }
0174 } else {
0175 external = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0176 internal = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0177 }
0178
0179 EgammaHcalIsolation thisHcalVar_ = EgammaHcalIsolation(external,
0180 outerCone_,
0181 internal,
0182 innerCone_,
0183 eThresHB_,
0184 etThresHB_,
0185 maxSeverityHB_,
0186 eThresHE_,
0187 etThresHE_,
0188 maxSeverityHE_,
0189 iEvent.get(hbheRecHitsTag_),
0190 iSetup.getData(caloGeometryToken_),
0191 iSetup.getData(hcalTopologyToken_),
0192 iSetup.getData(hcalChannelQualityToken_),
0193 iSetup.getData(hcalSevLvlComputerToken_),
0194 iSetup.getData(caloTowerConstituentsMapToken_));
0195
0196 if (useSingleTower_) {
0197 if (doEtSum_) {
0198 isol = thisHcalVar_.getHcalEtSumBc(recoEcalCandRef.get(), depth_);
0199 } else {
0200 isol = thisHcalVar_.getHcalESumBc(recoEcalCandRef.get(), depth_);
0201 }
0202 } else {
0203 if (doEtSum_) {
0204 isol = thisHcalVar_.getHcalEtSum(recoEcalCandRef.get(), depth_);
0205 } else {
0206 isol = thisHcalVar_.getHcalESum(recoEcalCandRef.get(), depth_);
0207 }
0208 }
0209
0210 if (doRhoCorrection_) {
0211 int iEA = -1;
0212 auto scEta = std::abs(recoEcalCandRef->superCluster()->eta());
0213 for (int bIt = absEtaLowEdges_.size() - 1; bIt > -1; bIt--) {
0214 if (scEta > absEtaLowEdges_.at(bIt)) {
0215 iEA = bIt;
0216 break;
0217 }
0218 }
0219 isol = isol - rho * effectiveAreas_.at(iEA);
0220 }
0221
0222 isoMap.insert(recoEcalCandRef, isol);
0223 }
0224
0225 iEvent.emplace(putToken_, isoMap);
0226 }
0227
0228 #include "FWCore/Framework/interface/MakerMacros.h"
0229 DEFINE_FWK_MODULE(EgammaHLTHcalVarProducerFromRecHit);