EgammaHLTProducers/plugins/EgammaHLTBcHcalIsolationProducersRegional.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaHLTProducers
0004 // Class:      EgammaHLTBcHcalIsolationProducersRegional
0005 //
0006 // Original Author:  Matteo Sani (UCSD)
0007 //         Created:  Thu Nov 24 11:38:00 CEST 2011
0008 //
0009
0010 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0011 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0012 #include "FWCore/Framework/interface/global/EDProducer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h"
0017 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
0018 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0019
0020 //this class produces either Hcal isolation or H for H/E  depending if doEtSum=true or false
0021 //H for H/E = towers behind SC, hcal isolation has these towers excluded
0022 //a rho correction can be applied
0023
0024 class EgammaHLTBcHcalIsolationProducersRegional : public edm::global::EDProducer<> {
0025 public:
0026   explicit EgammaHLTBcHcalIsolationProducersRegional(const edm::ParameterSet &);
0027
0028 public:
0029   void produce(edm::StreamID, edm::Event &, const edm::EventSetup &) const final;
0030   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0031
0032 private:
0033   const bool doEtSum_;
0034   const double etMin_;
0035   const double innerCone_;
0036   const double outerCone_;
0037   const int depth_;
0038   const bool useSingleTower_;
0039
0040   const bool doRhoCorrection_;
0041   const double rhoScale_;
0042   const double rhoMax_;
0043   const std::vector<double> effectiveAreas_;
0044   const std::vector<double> absEtaLowEdges_;
0045
0046   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateProducer_;
0047   const edm::EDGetTokenT<CaloTowerCollection> caloTowerProducer_;
0048   const edm::EDGetTokenT<double> rhoProducer_;
0049
0050   const edm::ESGetToken<CaloTowerConstituentsMap, CaloGeometryRecord> caloTowerConstituentsMapToken_;
0051
0052   const edm::EDPutTokenT<reco::RecoEcalCandidateIsolationMap> putToken_;
0053 };
0054
0055 EgammaHLTBcHcalIsolationProducersRegional::EgammaHLTBcHcalIsolationProducersRegional(const edm::ParameterSet &config)
0056     : doEtSum_(config.getParameter<bool>("doEtSum")),
0057       etMin_(config.getParameter<double>("etMin")),
0058       innerCone_(config.getParameter<double>("innerCone")),
0059       outerCone_(config.getParameter<double>("outerCone")),
0060       depth_(config.getParameter<int>("depth")),
0061       useSingleTower_(config.getParameter<bool>("useSingleTower")),
0062       doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0063       rhoScale_(config.getParameter<double>("rhoScale")),
0064       rhoMax_(config.getParameter<double>("rhoMax")),
0065       effectiveAreas_(config.getParameter<std::vector<double> >("effectiveAreas")),
0066       absEtaLowEdges_(config.getParameter<std::vector<double> >("absEtaLowEdges")),
0067       recoEcalCandidateProducer_(consumes(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
0068       caloTowerProducer_(consumes(config.getParameter<edm::InputTag>("caloTowerProducer"))),
0069       rhoProducer_(doRhoCorrection_ ? consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))
0070                                     : edm::EDGetTokenT<double>()),
0071       caloTowerConstituentsMapToken_{esConsumes()},
0072       putToken_{produces<reco::RecoEcalCandidateIsolationMap>()} {
0073   if (doRhoCorrection_) {
0074     if (absEtaLowEdges_.size() != effectiveAreas_.size()) {
0075       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0076     }
0077
0078     if (absEtaLowEdges_.at(0) != 0.0) {
0079       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0080     }
0081
0082     for (unsigned int aIt = 0; aIt < absEtaLowEdges_.size() - 1; aIt++) {
0083       if (!(absEtaLowEdges_.at(aIt) < absEtaLowEdges_.at(aIt + 1))) {
0084         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0085       }
0086     }
0087   }
0088 }
0089
0090 void EgammaHLTBcHcalIsolationProducersRegional::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0091   edm::ParameterSetDescription desc;
0092
0093   desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltRecoEcalCandidate"));
0094   desc.add<edm::InputTag>(("caloTowerProducer"), edm::InputTag("hltTowerMakerForAll"));
0095   desc.add<edm::InputTag>(("rhoProducer"), edm::InputTag("fixedGridRhoFastjetAllCalo"));
0096   desc.add<bool>(("doRhoCorrection"), false);
0097   desc.add<double>(("rhoMax"), 999999.);
0098   desc.add<double>(("rhoScale"), 1.0);
0099   desc.add<double>(("etMin"), -1.0);
0100   desc.add<double>(("innerCone"), 0);
0101   desc.add<double>(("outerCone"), 0.15);
0102   desc.add<int>(("depth"), -1);
0103   desc.add<bool>(("doEtSum"), false);
0104   desc.add<bool>(("useSingleTower"), false);
0105   desc.add<std::vector<double> >("effectiveAreas", {0.079, 0.25});  // 2016 post-ichep sinEle default
0106   desc.add<std::vector<double> >("absEtaLowEdges", {0.0, 1.479});   // Barrel, Endcap
0107   descriptions.add(("hltEgammaHLTBcHcalIsolationProducersRegional"), desc);
0108 }
0109
0110 void EgammaHLTBcHcalIsolationProducersRegional::produce(edm::StreamID,
0111                                                         edm::Event &iEvent,
0112                                                         const edm::EventSetup &iSetup) const {
0113   // Get the HLT filtered objects
0114   auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateProducer_);
0115
0116   double rho = 0.0;
0117
0118   if (doRhoCorrection_) {
0119     rho = iEvent.get(rhoProducer_);
0120     if (rho > rhoMax_) {
0121       rho = rhoMax_;
0122     }
0123     rho = rho * rhoScale_;
0124   }
0125
0126   auto const &caloTowers = iEvent.get(caloTowerProducer_);
0127   auto const &ctmaph = iSetup.getData(caloTowerConstituentsMapToken_);
0128
0129   const EgammaTowerIsolation towerIso1(outerCone_, 0., etMin_, 1, &caloTowers);
0130   const EgammaTowerIsolation towerIso2(outerCone_, 0., etMin_, 2, &caloTowers);
0131
0132   reco::RecoEcalCandidateIsolationMap isoMap(recoEcalCandHandle);
0133
0134   for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoEcalCandHandle->size(); iRecoEcalCand++) {
0135     reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand);
0136
0137     float isol = 0;
0138
0139     auto towersBehindCluster =
0140         useSingleTower_ ? egamma::towersOf(*(recoEcalCandRef->superCluster()), ctmaph) : std::vector<CaloTowerDetId>{};
0141
0142     if (doEtSum_) {  //calculate hcal isolation excluding the towers behind the cluster which will be used for H for H/E
0143       const EgammaTowerIsolation isolAlgo(outerCone_, innerCone_, etMin_, depth_, &caloTowers);
0144       if (useSingleTower_) {
0145         // towersBehindCluster are excluded from the isolation sum
0146         isol = isolAlgo.getTowerEtSum(recoEcalCandRef.get(), &towersBehindCluster);
0147       } else {
0148         isol = isolAlgo.getTowerEtSum(recoEcalCandRef.get());
0149       }
0150
0151     } else {  //calcuate H for H/E
0152       if (useSingleTower_)
0153         isol = egamma::depth1HcalESum(towersBehindCluster, caloTowers) +
0154                egamma::depth2HcalESum(towersBehindCluster, caloTowers);
0155       else {
0156         auto const &sc = recoEcalCandRef->superCluster().get();
0157         isol = towerIso1.getTowerESum(sc) + towerIso2.getTowerESum(sc);
0158       }
0159     }
0160
0161     if (doRhoCorrection_) {
0162       int iEA = -1;
0163       auto scEta = std::abs(recoEcalCandRef->superCluster()->eta());
0164       for (int bIt = absEtaLowEdges_.size() - 1; bIt > -1; bIt--) {
0165         if (scEta >= absEtaLowEdges_[bIt]) {
0166           iEA = bIt;
0167           break;
0168         }
0169       }
0170       isol = isol - rho * effectiveAreas_[iEA];
0171     }
0172     isoMap.insert(recoEcalCandRef, isol);
0173   }
0174
0175   iEvent.emplace(putToken_, isoMap);
0176 }
0177
0178 #include "FWCore/Framework/interface/MakerMacros.h"
0179 DEFINE_FWK_MODULE(EgammaHLTBcHcalIsolationProducersRegional);