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
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   //eThresHB/HE are from RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py
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});  // 2016 post-ichep sinEle default
0138   desc.add<std::vector<double> >("absEtaLowEdges", {0.0, 1.479});   // Barrel, Endcap
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_) {  //this is single tower based H/E
0168         external = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0169         internal = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0170       } else {  //this is cone-based HCAL isolation with single tower based footprint removal
0171         external = EgammaHcalIsolation::InclusionRule::withinConeAroundCluster;
0172         internal = EgammaHcalIsolation::InclusionRule::isBehindClusterSeed;
0173       }
0174     } else {  //useSingleTower_=False means H/E is cone-based
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_) {  //this is cone-based HCAL isolation with single tower based footprint removal
0198         isol = thisHcalVar_.getHcalEtSumBc(recoEcalCandRef.get(), depth_);  //depth=0 means all depths
0199       } else {                                                              //this is single tower based H/E
0200         isol = thisHcalVar_.getHcalESumBc(recoEcalCandRef.get(), depth_);   //depth=0 means all depths
0201       }
0202     } else {           //useSingleTower_=False means H/E is cone-based.
0203       if (doEtSum_) {  //hcal iso
0204         isol = thisHcalVar_.getHcalEtSum(recoEcalCandRef.get(), depth_);  //depth=0 means all depths
0205       } else {  // doEtSum_=False means sum up energy, this is for H/E
0206         isol = thisHcalVar_.getHcalESum(recoEcalCandRef.get(), depth_);  //depth=0 means all depths
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);