File indexing completed on 2023-03-17 11:17:41
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
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
0039 explicit EgammaIsoDetIdCollectionProducer(const edm::ParameterSet&);
0040
0041 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0042
0043 private:
0044
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
0093 produces<DetIdCollection>(interestingDetIdCollection_);
0094 }
0095
0096
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
0123 auto detIdCollection = std::make_unique<DetIdCollection>();
0124
0125 if (doubleConeSel_) {
0126 for (auto const& emObj : emObjects) {
0127
0128 if (emObj.et() < etCandCut_)
0129 continue;
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;
0135
0136 double et =
0137 recIt.energy() * caloGeom->getPosition(recIt.detid()).perp() / caloGeom->getPosition(recIt.detid()).mag();
0138
0139 if (et < etCut_)
0140 return;
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
0159 if (!recIt.checkFlag(EcalRecHit::kGood)) {
0160 if (recIt.checkFlags(flagsexclEB_)) {
0161 return;
0162 }
0163 }
0164 } else {
0165
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 });
0176
0177 }
0178
0179 }
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);