EgammaIsolationAlgos/plugins/EgammaIsoDetIdCollectionProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaIsoDetIdCollectionProducer
0004 // Class:      EgammaIsoDetIdCollectionProducer
0005 //
0006 /**\class EgammaIsoDetIdCollectionProducer
0007 Original author: Matthew LeBourgeois PH/CMG
0008 Modified from :
0009 RecoEcal/EgammaClusterProducers/{src,interface}/InterestingDetIdCollectionProducer.{h,cc}
0010 by Paolo Meridiani PH/CMG
0011
0012 Implementation:
0013  <Notes on implementation>
0014 */
0015
0016 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0017 #include "DataFormats/DetId/interface/DetIdCollection.h"
0018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0019 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0020 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0021 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0022 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030 #include "RecoCaloTools/Selectors/interface/CaloDualConeSelector.h"
0031 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0032 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0033
0034 template <class T1>
0035 class EgammaIsoDetIdCollectionProducer : public edm::global::EDProducer<> {
0036 public:
0037   typedef std::vector<T1> T1Collection;
0038   //! ctor
0039   explicit EgammaIsoDetIdCollectionProducer(const edm::ParameterSet&);
0040   //! producer
0041   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0042
0043 private:
0044   // ----------member data ---------------------------
0045   edm::EDGetTokenT<EcalRecHitCollection> recHitsToken_;
0046   edm::EDGetTokenT<T1Collection> emObjectToken_;
0047   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0048   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLvToken_;
0049   edm::InputTag recHitsLabel_;
0050   edm::InputTag emObjectLabel_;
0051   double energyCut_;
0052   double etCut_;
0053   double etCandCut_;
0054   double outerRadius_;
0055   double innerRadius_;
0056   std::string interestingDetIdCollection_;
0057
0058   std::vector<int> severitiesexclEB_;
0059   std::vector<int> severitiesexclEE_;
0060   std::vector<int> flagsexclEB_;
0061   std::vector<int> flagsexclEE_;
0062 };
0063
0064 template <class T1>
0065 EgammaIsoDetIdCollectionProducer<T1>::EgammaIsoDetIdCollectionProducer(const edm::ParameterSet& iConfig)
0066     : recHitsToken_{consumes(iConfig.getParameter<edm::InputTag>("recHitsLabel"))},
0067       emObjectToken_{consumes(iConfig.getParameter<edm::InputTag>("emObjectLabel"))},
0068       caloGeometryToken_(esConsumes()),
0069       sevLvToken_(esConsumes()),
0070       recHitsLabel_(iConfig.getParameter<edm::InputTag>("recHitsLabel")),
0071       emObjectLabel_(iConfig.getParameter<edm::InputTag>("emObjectLabel")),
0072       energyCut_(iConfig.getParameter<double>("energyCut")),
0073       etCut_(iConfig.getParameter<double>("etCut")),
0074       etCandCut_(iConfig.getParameter<double>("etCandCut")),
0075       outerRadius_(iConfig.getParameter<double>("outerRadius")),
0076       innerRadius_(iConfig.getParameter<double>("innerRadius")),
0077       interestingDetIdCollection_(iConfig.getParameter<std::string>("interestingDetIdCollection")) {
0078   auto const& flagnamesEB = iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB");
0079   auto const& flagnamesEE = iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE");
0080
0081   flagsexclEB_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
0082   flagsexclEE_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
0083
0084   auto const& severitynamesEB = iConfig.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEB");
0085
0086   severitiesexclEB_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
0087
0088   auto const& severitynamesEE = iConfig.getParameter<std::vector<std::string>>("RecHitSeverityToBeExcludedEE");
0089
0090   severitiesexclEE_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
0091
0092   //register your products
0093   produces<DetIdCollection>(interestingDetIdCollection_);
0094 }
0095
0096 // ------------ method called to produce the data  ------------
0097 template <class T1>
0098 void EgammaIsoDetIdCollectionProducer<T1>::produce(edm::StreamID,
0099                                                    edm::Event& iEvent,
0100                                                    const edm::EventSetup& iSetup) const {
0101   using namespace edm;
0102   using namespace std;
0103
0104   auto const& emObjects = iEvent.get(emObjectToken_);
0105   auto const& ecalRecHits = iEvent.get(recHitsToken_);
0106
0107   edm::ESHandle<CaloGeometry> pG = iSetup.getHandle(caloGeometryToken_);
0108   const CaloGeometry* caloGeom = pG.product();
0109
0110   edm::ESHandle<EcalSeverityLevelAlgo> sevlv = iSetup.getHandle(sevLvToken_);
0111   const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
0112
0113   std::unique_ptr<CaloDualConeSelector<EcalRecHit>> doubleConeSel_ = nullptr;
0114   if (recHitsLabel_.instance() == "EcalRecHitsEB") {
0115     doubleConeSel_ =
0116         std::make_unique<CaloDualConeSelector<EcalRecHit>>(innerRadius_, outerRadius_, &*pG, DetId::Ecal, EcalBarrel);
0117   } else if (recHitsLabel_.instance() == "EcalRecHitsEE") {
0118     doubleConeSel_ =
0119         std::make_unique<CaloDualConeSelector<EcalRecHit>>(innerRadius_, outerRadius_, &*pG, DetId::Ecal, EcalEndcap);
0120   }
0121
0122   //Create empty output collections
0123   auto detIdCollection = std::make_unique<DetIdCollection>();
0124
0125   if (doubleConeSel_) {                    //if cone selector was created
0126     for (auto const& emObj : emObjects) {  //Loop over candidates
0127
0128       if (emObj.et() < etCandCut_)
0129         continue;  //don't calculate if object hasn't enough energy
0130
0131       GlobalPoint pclu(emObj.caloPosition().x(), emObj.caloPosition().y(), emObj.caloPosition().z());
0132       doubleConeSel_->selectCallback(pclu, ecalRecHits, [&](const EcalRecHit& recIt) {
0133         if (recIt.energy() < energyCut_)
0134           return;  //dont fill if below E noise value
0135
0136         double et =
0137             recIt.energy() * caloGeom->getPosition(recIt.detid()).perp() / caloGeom->getPosition(recIt.detid()).mag();
0138
0139         if (et < etCut_)
0140           return;  //dont fill if below ET noise value
0141
0142         bool isBarrel = false;
0143         if (fabs(caloGeom->getPosition(recIt.detid()).eta() < 1.479))
0144           isBarrel = true;
0145
0146         int severityFlag = sevLevel->severityLevel(recIt.detid(), ecalRecHits);
0147         if (isBarrel) {
0148           auto sit = std::find(severitiesexclEB_.begin(), severitiesexclEB_.end(), severityFlag);
0149           if (sit != severitiesexclEB_.end())
0150             return;
0151         } else {
0152           auto sit = std::find(severitiesexclEE_.begin(), severitiesexclEE_.end(), severityFlag);
0153           if (sit != severitiesexclEE_.end())
0154             return;
0155         }
0156
0157         if (isBarrel) {
0158           // new rechit flag checks
0159           if (!recIt.checkFlag(EcalRecHit::kGood)) {
0160             if (recIt.checkFlags(flagsexclEB_)) {
0161               return;
0162             }
0163           }
0164         } else {
0165           // new rechit flag checks
0166           if (!recIt.checkFlag(EcalRecHit::kGood)) {
0167             if (recIt.checkFlags(flagsexclEE_)) {
0168               return;
0169             }
0170           }
0171         }
0172
0173         if (std::find(detIdCollection->begin(), detIdCollection->end(), recIt.detid()) == detIdCollection->end())
0174           detIdCollection->push_back(recIt.detid());
0175       });  //end rechits
0176
0177     }  //end candidates
0178
0179   }  //end if cone selector was created
0180
0181   iEvent.put(std::move(detIdCollection), interestingDetIdCollection_);
0182 }
0183
0184 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0185 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0186 #include "FWCore/Framework/interface/MakerMacros.h"
0187
0188 using EleIsoDetIdCollectionProducer = EgammaIsoDetIdCollectionProducer<reco::GsfElectron>;
0189 using GamIsoDetIdCollectionProducer = EgammaIsoDetIdCollectionProducer<reco::Photon>;
0190
0191 DEFINE_FWK_MODULE(GamIsoDetIdCollectionProducer);
0192 DEFINE_FWK_MODULE(EleIsoDetIdCollectionProducer);