EgammaHLTProducers/plugins/EgammaHLTHcalVarProducerFromRecHit.cc

0001 // Class: EgammaHLTHcalVarProducerFromRecHit
0002
0003 /*
0004
0005 Author: Swagata Mukherjee
0006
0007 Date: August 2021
0008
0009 This class is similar to the existing class EgammaHLTBcHcalIsolationProducersRegional,
0010 but the new feature in this code is that the HCAL recHits are used instead of the
0011 calotowers which is expected to be phased out sometime in Run3.
0012 The old class can also be used until calotowers stay. After that, one need to switch to this new one.
0013
0014 As the old producer code, this one also produces either Hcal isolation or H for H/E depending if doEtSum=true or false.
0015 H for H/E = either a single HCAL tower behind SC, or towers in a cone, and hcal isolation has these tower(s) excluded.
0016 A rho correction can be applied
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 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0036 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0037 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0038
0039 class EgammaHLTHcalVarProducerFromRecHit : public edm::global::EDProducer<> {
0040 public:
0041   explicit EgammaHLTHcalVarProducerFromRecHit(const edm::ParameterSet &);
0042
0043 public:
0044   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const final;
0045   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0046
0047 private:
0048   const bool doEtSum_;
0049   const EgammaHcalIsolation::arrayHB eThresHB_;
0050   const EgammaHcalIsolation::arrayHB etThresHB_;
0051   const EgammaHcalIsolation::arrayHE eThresHE_;
0052   const EgammaHcalIsolation::arrayHE etThresHE_;
0053   const double innerCone_;
0054   const double outerCone_;
0055   const int depth_;
0056   const int maxSeverityHB_;
0057   const int maxSeverityHE_;
0058   const bool useSingleTower_;
0059   const bool doRhoCorrection_;
0060   const double rhoScale_;
0061   const double rhoMax_;
0062   const std::vector<double> effectiveAreas_;
0063   const std::vector<double> absEtaLowEdges_;
0064   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateProducer_;
0065   const edm::EDGetTokenT<HBHERecHitCollection> hbheRecHitsTag_;
0066   const edm::EDGetTokenT<double> rhoProducer_;
0067   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0068   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalTopologyToken_;
0069   const edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> hcalChannelQualityToken_;
0070   const edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> hcalSevLvlComputerToken_;
0071   const edm::ESGetToken<CaloTowerConstituentsMap, CaloGeometryRecord> caloTowerConstituentsMapToken_;
0072   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> putToken_;
0073
0074   //Get HCAL thresholds from GT
0075   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0076   bool cutsFromDB;
0077 };
0078
0079 EgammaHLTHcalVarProducerFromRecHit::EgammaHLTHcalVarProducerFromRecHit(const edm::ParameterSet &config)
0080     : doEtSum_(config.getParameter<bool>("doEtSum")),
0081       eThresHB_(config.getParameter<EgammaHcalIsolation::arrayHB>("eThresHB")),
0082       etThresHB_(config.getParameter<EgammaHcalIsolation::arrayHB>("etThresHB")),
0083       eThresHE_(config.getParameter<EgammaHcalIsolation::arrayHE>("eThresHE")),
0084       etThresHE_(config.getParameter<EgammaHcalIsolation::arrayHE>("etThresHE")),
0085       innerCone_(config.getParameter<double>("innerCone")),
0086       outerCone_(config.getParameter<double>("outerCone")),
0087       depth_(config.getParameter<int>("depth")),
0088       maxSeverityHB_(config.getParameter<int>("maxSeverityHB")),
0089       maxSeverityHE_(config.getParameter<int>("maxSeverityHE")),
0090       useSingleTower_(config.getParameter<bool>("useSingleTower")),
0091       doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0092       rhoScale_(config.getParameter<double>("rhoScale")),
0093       rhoMax_(config.getParameter<double>("rhoMax")),
0094       effectiveAreas_(config.getParameter<std::vector<double> >("effectiveAreas")),
0095       absEtaLowEdges_(config.getParameter<std::vector<double> >("absEtaLowEdges")),
0096       recoEcalCandidateProducer_(consumes(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0097       hbheRecHitsTag_(consumes(config.getParameter<edm::InputTag>("hbheRecHitsTag"))),
0098       rhoProducer_(doRhoCorrection_ ? consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))
0099                                     : edm::EDGetTokenT<double>()),
0100       caloGeometryToken_{esConsumes()},
0101       hcalTopologyToken_{esConsumes()},
0102       hcalChannelQualityToken_{esConsumes(edm::ESInputTag("", "withTopo"))},
0103       hcalSevLvlComputerToken_{esConsumes()},
0104       caloTowerConstituentsMapToken_{esConsumes()},
0105       putToken_{produces<reco::RecoEcalCandidateIsolationMap>()},
0106       cutsFromDB(
0107           config.getParameter<bool>("usePFThresholdsFromDB")) {  //Retrieve HCAL PF thresholds - from config or from DB
0108   if (doRhoCorrection_) {
0109     if (absEtaLowEdges_.size() != effectiveAreas_.size()) {
0110       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0111     }
0112
0113     if (absEtaLowEdges_.at(0) != 0.0) {
0114       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0115     }
0116
0117     for (unsigned int aIt = 0; aIt < absEtaLowEdges_.size() - 1; aIt++) {
0118       if (!(absEtaLowEdges_.at(aIt) < absEtaLowEdges_.at(aIt + 1))) {
0119         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0120       }
0121     }
0122   }
0123
0124   if (cutsFromDB) {
0125     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
0126   }
0127 }
0128
0129 void EgammaHLTHcalVarProducerFromRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0130   edm::ParameterSetDescription desc;
0131
0132   desc.add<edm::InputTag>("recoEcalCandidateProducer", edm::InputTag("hltRecoEcalCandidate"));
0133   desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0134   desc.add<edm::InputTag>("hbheRecHitsTag", edm::InputTag("hltHbhereco"));
0135   desc.add<bool>("doRhoCorrection", false);
0136   desc.add<double>("rhoMax", 999999.);
0137   desc.add<double>(("rhoScale"), 1.0);
0138   //eThresHB/HE are from RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py
0139   desc.add<std::vector<double> >("eThresHB", {0.1, 0.2, 0.3, 0.3});
0140   desc.add<std::vector<double> >("etThresHB", {0, 0, 0, 0});
0141   desc.add<std::vector<double> >("eThresHE", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0142   desc.add<std::vector<double> >("etThresHE", {0, 0, 0, 0, 0, 0, 0});
0143   desc.add<bool>("usePFThresholdsFromDB", true);
0144   desc.add<double>("innerCone", 0);
0145   desc.add<double>("outerCone", 0.14);
0146   desc.add<int>("depth", 0);
0147   desc.add<int>("maxSeverityHB", 9);
0148   desc.add<int>("maxSeverityHE", 9);
0149   desc.add<bool>("doEtSum", false);
0150   desc.add<bool>("useSingleTower", false);
0151   desc.add<std::vector<double> >("effectiveAreas", {0.079, 0.25});  // 2016 post-ichep sinEle default
0152   desc.add<std::vector<double> >("absEtaLowEdges", {0.0, 1.479});   // Barrel, Endcap
0153   descriptions.add("hltEgammaHLTHcalVarProducerFromRecHit", desc);
0154 }
0155
0156 void EgammaHLTHcalVarProducerFromRecHit::produce(edm::StreamID,
0157                                                  edm::Event &iEvent,
0158                                                  const edm::EventSetup &iSetup) const {
0159   auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateProducer_);
0160
0161   double rho = 0.0;
0162
0163   if (doRhoCorrection_) {
0164     rho = iEvent.get(rhoProducer_);
0165     if (rho > rhoMax_) {
0166       rho = rhoMax_;
0167     }
0168     rho = rho * rhoScale_;
0169   }
0170
0171   reco::RecoEcalCandidateIsolationMap isoMap(recoEcalCandHandle);
0172
0173   for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); iRecoEcalCand++) {
0174     reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
0175
0176     float isol = 0;
0177     EgammaHcalIsolation::InclusionRule external;
0178     EgammaHcalIsolation::InclusionRule internal;
0179
0180     if (useSingleTower_) {
0181       if (!doEtSum_) {  //this is single tower based H/E
0182         external = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0183         internal = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0184       } else {  //this is cone-based HCAL isolation with single tower based footprint removal
0185         external = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0186         internal = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0187       }
0188     } else {  //useSingleTower_=False means H/E is cone-based
0189       external = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0190       internal = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0191     }
0192
0193     EgammaHcalIsolation thisHcalVar_ = EgammaHcalIsolation(external,
0194                                                            outerCone_,
0195                                                            internal,
0196                                                            innerCone_,
0197                                                            eThresHB_,
0198                                                            etThresHB_,
0199                                                            maxSeverityHB_,
0200                                                            eThresHE_,
0201                                                            etThresHE_,
0202                                                            maxSeverityHE_,
0203                                                            iEvent.get(hbheRecHitsTag_),
0204                                                            iSetup.getData(caloGeometryToken_),
0205                                                            iSetup.getData(hcalTopologyToken_),
0206                                                            iSetup.getData(hcalChannelQualityToken_),
0207                                                            iSetup.getData(hcalSevLvlComputerToken_),
0208                                                            iSetup.getData(caloTowerConstituentsMapToken_));
0209     const HcalPFCuts *hcalCuts{nullptr};
0210     if (cutsFromDB) {
0211       hcalCuts = &iSetup.getData(hcalCutsToken_);
0212     }
0213
0214     if (useSingleTower_) {
0215       if (doEtSum_) {  //this is cone-based HCAL isolation with single tower based footprint removal
0216         isol = thisHcalVar_.getHcalEtSumBc(recoEcalCandRef.get(), depth_, hcalCuts);  //depth=0 means all depths
0217       } else {                                                                        //this is single tower based H/E
0218         isol = thisHcalVar_.getHcalESumBc(recoEcalCandRef.get(), depth_, hcalCuts);   //depth=0 means all depths
0219       }
0220     } else {           //useSingleTower_=False means H/E is cone-based.
0221       if (doEtSum_) {  //hcal iso
0222         isol = thisHcalVar_.getHcalEtSum(recoEcalCandRef.get(), depth_, hcalCuts);  //depth=0 means all depths
0223       } else {  // doEtSum_=False means sum up energy, this is for H/E
0224         isol = thisHcalVar_.getHcalESum(recoEcalCandRef.get(), depth_, hcalCuts);  //depth=0 means all depths
0225       }
0226     }
0227
0228     if (doRhoCorrection_) {
0229       int iEA = -1;
0230       auto scEta = std::abs(recoEcalCandRef->superCluster()->eta());
0231       for (int bIt = absEtaLowEdges_.size() - 1; bIt > -1; bIt--) {
0232         if (scEta >= absEtaLowEdges_[bIt]) {
0233           iEA = bIt;
0234           break;
0235         }
0236       }
0237       isol = isol - rho * effectiveAreas_[iEA];
0238     }
0239
0240     isoMap.insert(recoEcalCandRef, isol);
0241   }
0242
0243   iEvent.emplace(putToken_, isoMap);
0244 }
0245
0246 #include "FWCore/Framework/interface/MakerMacros.h"
0247 DEFINE_FWK_MODULE(EgammaHLTHcalVarProducerFromRecHit);