File indexing completed on 2024-04-06 12:25:07
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/ESGetToken.h"
0010
0011 #include "DataFormats/Common/interface/ValueMap.h"
0012 #include "DataFormats/Common/interface/View.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014
0015 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0016 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0017 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0018 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0019 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0020
0021 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0022 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0023
0024 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
0026
0027 #include <memory>
0028 #include <vector>
0029
0030 using IsolationCalculators = std::vector<std::unique_ptr<EleTkIsolFromCands>>;
0031
0032
0033 class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> {
0034 private:
0035
0036 template <typename T>
0037 struct DualToken {
0038 edm::EDGetTokenT<T> aod;
0039 edm::EDGetTokenT<T> miniAOD;
0040 };
0041 class DataFormat {
0042 public:
0043 enum Format { AUTO = 0, AOD = 1, MINIAOD = 2 };
0044
0045 private:
0046 int data_;
0047
0048 public:
0049 DataFormat(int val) : data_(val) {}
0050 bool tryAOD() const { return data_ == AUTO || data_ == AOD; }
0051 bool tryMiniAOD() const { return data_ == AUTO || data_ == MINIAOD; }
0052 int operator()() const { return data_; }
0053 };
0054
0055 public:
0056 explicit ElectronHEEPIDValueMapProducer(const edm::ParameterSet&);
0057 ~ElectronHEEPIDValueMapProducer() override;
0058
0059 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062 void produce(edm::Event&, const edm::EventSetup&) override;
0063
0064 template <typename T>
0065 static void writeValueMap(edm::Event& iEvent,
0066 const edm::Handle<edm::View<reco::GsfElectron>>& handle,
0067 const std::vector<T>& values,
0068 const std::string& label);
0069
0070 static int nrSaturatedCrysIn5x5(const reco::GsfElectron& ele,
0071 edm::Handle<EcalRecHitCollection>& ebHits,
0072 edm::Handle<EcalRecHitCollection>& eeHits,
0073 edm::ESHandle<CaloTopology>& caloTopo);
0074
0075 static float calTrkIso(const reco::GsfElectron& ele, IsolationCalculators const& isolationCalculators);
0076
0077 template <typename T>
0078 void setToken(edm::EDGetTokenT<T>& token, edm::InputTag tag) {
0079 token = consumes<T>(tag);
0080 }
0081 template <typename T>
0082 void setToken(edm::EDGetTokenT<T>& token, const edm::ParameterSet& iPara, const std::string& tag) {
0083 token = consumes<T>(iPara.getParameter<edm::InputTag>(tag));
0084 }
0085 template <typename T>
0086 void setToken(std::vector<edm::EDGetTokenT<T>>& tokens, const edm::ParameterSet& iPara, const std::string& tagName) {
0087 auto tags = iPara.getParameter<std::vector<edm::InputTag>>(tagName);
0088 for (auto& tag : tags) {
0089 edm::EDGetTokenT<T> token;
0090 setToken(token, tag);
0091 tokens.push_back(token);
0092 }
0093 }
0094 template <typename T>
0095 void setToken(DualToken<T>& token,
0096 const edm::ParameterSet& iPara,
0097 const std::string& tagAOD,
0098 const std::string& tagMiniAOD,
0099 DataFormat format) {
0100 if (format.tryAOD())
0101 token.aod = consumes<T>(iPara.getParameter<edm::InputTag>(tagAOD));
0102 if (format.tryMiniAOD())
0103 token.miniAOD = consumes<T>(iPara.getParameter<edm::InputTag>(tagMiniAOD));
0104 }
0105 template <typename T>
0106 void setToken(std::vector<DualToken<T>>& tokens,
0107 const edm::ParameterSet& iPara,
0108 const std::string& tagAOD,
0109 const std::string& tagMiniAOD,
0110 DataFormat format) {
0111 auto tagsAOD = iPara.getParameter<std::vector<edm::InputTag>>(tagAOD);
0112 auto tagsMiniAOD = iPara.getParameter<std::vector<edm::InputTag>>(tagMiniAOD);
0113 size_t maxSize = std::max(tagsAOD.size(), tagsMiniAOD.size());
0114 tokens.clear();
0115 tokens.resize(maxSize);
0116 if (format.tryAOD()) {
0117 for (size_t tagNr = 0; tagNr < tagsAOD.size(); tagNr++) {
0118 setToken(tokens[tagNr].aod, tagsAOD[tagNr]);
0119 }
0120 }
0121 if (format.tryMiniAOD()) {
0122 for (size_t tagNr = 0; tagNr < tagsMiniAOD.size(); tagNr++) {
0123 setToken(tokens[tagNr].miniAOD, tagsMiniAOD[tagNr]);
0124 }
0125 }
0126 }
0127
0128 template <typename T>
0129 static edm::Handle<T> getHandle(const edm::Event& iEvent, const DualToken<T>& token) {
0130 edm::Handle<T> handle;
0131 if (!token.aod.isUninitialized())
0132 iEvent.getByToken(token.aod, handle);
0133 if (!handle.isValid() && !token.miniAOD.isUninitialized())
0134 iEvent.getByToken(token.miniAOD, handle);
0135 return handle;
0136 }
0137
0138 template <typename T>
0139 static std::vector<edm::Handle<T>> getHandles(const edm::Event& iEvent, const std::vector<DualToken<T>>& tokens) {
0140 std::vector<edm::Handle<T>> handles(tokens.size());
0141 if (tokens.empty())
0142 return handles;
0143 if (!tokens[0].aod.isUninitialized())
0144 iEvent.getByToken(tokens[0].aod, handles[0]);
0145 bool isAOD = handles[0].isValid();
0146 if (!isAOD && !tokens[0].miniAOD.isUninitialized())
0147 iEvent.getByToken(tokens[0].miniAOD, handles[0]);
0148
0149 for (size_t tokenNr = 1; tokenNr < tokens.size(); tokenNr++) {
0150 auto token = isAOD ? tokens[tokenNr].aod : tokens[tokenNr].miniAOD;
0151 if (!token.isUninitialized())
0152 iEvent.getByToken(token, handles[tokenNr]);
0153 }
0154 return handles;
0155 }
0156
0157 template <typename T>
0158 static bool isEventAOD(const edm::Event& iEvent, const DualToken<T>& token) {
0159 edm::Handle<T> handle;
0160 if (!token.aod.isUninitialized())
0161 iEvent.getByToken(token.aod, handle);
0162 if (handle.isValid())
0163 return true;
0164 else
0165 return false;
0166 }
0167
0168 DualToken<EcalRecHitCollection> ebRecHitToken_;
0169 DualToken<EcalRecHitCollection> eeRecHitToken_;
0170 DualToken<edm::View<reco::GsfElectron>> eleToken_;
0171 std::vector<DualToken<pat::PackedCandidateCollection>> candTokens_;
0172 edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0173 edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopoToken_;
0174
0175 EleTkIsolFromCands::Configuration trkIsoCalcCfg_;
0176 EleTkIsolFromCands::Configuration trkIso04CalcCfg_;
0177 bool makeTrkIso04_;
0178 DataFormat dataFormat_;
0179 std::vector<EleTkIsolFromCands::PIDVeto> candVetosAOD_;
0180 std::vector<EleTkIsolFromCands::PIDVeto> candVetosMiniAOD_;
0181
0182 static const std::string eleTrkPtIsoLabel_;
0183 static const std::string eleTrkPtIso04Label_;
0184 static const std::string eleNrSaturateIn5x5Label_;
0185 };
0186
0187 const std::string ElectronHEEPIDValueMapProducer::eleTrkPtIsoLabel_ = "eleTrkPtIso";
0188 const std::string ElectronHEEPIDValueMapProducer::eleTrkPtIso04Label_ = "eleTrkPtIso04";
0189 const std::string ElectronHEEPIDValueMapProducer::eleNrSaturateIn5x5Label_ = "eleNrSaturateIn5x5";
0190
0191 ElectronHEEPIDValueMapProducer::ElectronHEEPIDValueMapProducer(const edm::ParameterSet& iConfig)
0192 : trkIsoCalcCfg_(iConfig.getParameter<edm::ParameterSet>("trkIsoConfig")),
0193 trkIso04CalcCfg_(iConfig.getParameter<edm::ParameterSet>("trkIso04Config")),
0194 makeTrkIso04_(iConfig.getParameter<bool>("makeTrkIso04")),
0195 dataFormat_(iConfig.getParameter<int>("dataFormat")) {
0196 setToken(ebRecHitToken_, iConfig, "ebRecHitsAOD", "ebRecHitsMiniAOD", dataFormat_);
0197 setToken(eeRecHitToken_, iConfig, "eeRecHitsAOD", "eeRecHitsMiniAOD", dataFormat_);
0198 setToken(eleToken_, iConfig, "elesAOD", "elesMiniAOD", dataFormat_);
0199 setToken(candTokens_, iConfig, "candsAOD", "candsMiniAOD", dataFormat_);
0200 setToken(beamSpotToken_, iConfig, "beamSpot");
0201 caloTopoToken_ = esConsumes();
0202
0203 auto fillVetos = [](const auto& in, auto& out) {
0204 std::transform(in.begin(), in.end(), std::back_inserter(out), EleTkIsolFromCands::pidVetoFromStr);
0205 };
0206
0207 fillVetos(iConfig.getParameter<std::vector<std::string>>("candVetosAOD"), candVetosAOD_);
0208 if (candVetosAOD_.size() != iConfig.getParameter<std::vector<edm::InputTag>>("candsAOD").size()) {
0209 throw cms::Exception("ConfigError") << " Error candVetosAOD should be the same size as candsAOD " << std::endl;
0210 }
0211
0212 fillVetos(iConfig.getParameter<std::vector<std::string>>("candVetosMiniAOD"), candVetosMiniAOD_);
0213 if (candVetosMiniAOD_.size() != iConfig.getParameter<std::vector<edm::InputTag>>("candsMiniAOD").size()) {
0214 throw cms::Exception("ConfigError") << " Error candVetosMiniAOD should be the same size as candsMiniAOD "
0215 << std::endl;
0216 }
0217
0218 produces<edm::ValueMap<float>>(eleTrkPtIsoLabel_);
0219 if (makeTrkIso04_)
0220 produces<edm::ValueMap<float>>(eleTrkPtIso04Label_);
0221 produces<edm::ValueMap<int>>(eleNrSaturateIn5x5Label_);
0222 }
0223
0224 ElectronHEEPIDValueMapProducer::~ElectronHEEPIDValueMapProducer() {}
0225
0226 void ElectronHEEPIDValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0227 auto eleHandle = getHandle(iEvent, eleToken_);
0228 auto ebRecHitHandle = getHandle(iEvent, ebRecHitToken_);
0229 auto eeRecHitHandle = getHandle(iEvent, eeRecHitToken_);
0230 auto beamSpotHandle = iEvent.getHandle(beamSpotToken_);
0231
0232 bool isAOD = isEventAOD(iEvent, eleToken_);
0233 const auto& candVetos = isAOD ? candVetosAOD_ : candVetosMiniAOD_;
0234
0235 IsolationCalculators trkIsoCalcs;
0236 IsolationCalculators trkIso04Calcs;
0237
0238 {
0239 int handleNr = 0;
0240 for (auto const& handle : getHandles(iEvent, candTokens_)) {
0241 if (handle.isValid()) {
0242 trkIsoCalcs.emplace_back(std::make_unique<EleTkIsolFromCands>(trkIsoCalcCfg_, *handle, candVetos[handleNr]));
0243 if (makeTrkIso04_) {
0244 trkIso04Calcs.emplace_back(
0245 std::make_unique<EleTkIsolFromCands>(trkIso04CalcCfg_, *handle, candVetos[handleNr]));
0246 }
0247 }
0248 ++handleNr;
0249 }
0250 }
0251
0252 edm::ESHandle<CaloTopology> caloTopoHandle = iSetup.getHandle(caloTopoToken_);
0253
0254 std::vector<float> eleTrkPtIso;
0255 std::vector<float> eleTrkPtIso04;
0256 std::vector<int> eleNrSaturateIn5x5;
0257 for (auto const& ele : *eleHandle) {
0258 eleTrkPtIso.push_back(calTrkIso(ele, trkIsoCalcs));
0259 if (makeTrkIso04_) {
0260 eleTrkPtIso04.push_back(calTrkIso(ele, trkIso04Calcs));
0261 }
0262 eleNrSaturateIn5x5.push_back(nrSaturatedCrysIn5x5(ele, ebRecHitHandle, eeRecHitHandle, caloTopoHandle));
0263 }
0264
0265 writeValueMap(iEvent, eleHandle, eleTrkPtIso, eleTrkPtIsoLabel_);
0266 if (makeTrkIso04_)
0267 writeValueMap(iEvent, eleHandle, eleTrkPtIso04, eleTrkPtIso04Label_);
0268 writeValueMap(iEvent, eleHandle, eleNrSaturateIn5x5, eleNrSaturateIn5x5Label_);
0269 }
0270
0271 int ElectronHEEPIDValueMapProducer::nrSaturatedCrysIn5x5(const reco::GsfElectron& ele,
0272 edm::Handle<EcalRecHitCollection>& ebHits,
0273 edm::Handle<EcalRecHitCollection>& eeHits,
0274 edm::ESHandle<CaloTopology>& caloTopo) {
0275 DetId id = ele.superCluster()->seed()->seed();
0276 auto recHits = id.subdetId() == EcalBarrel ? ebHits.product() : eeHits.product();
0277 return noZS::EcalClusterTools::nrSaturatedCrysIn5x5(id, recHits, caloTopo.product());
0278 }
0279
0280 float ElectronHEEPIDValueMapProducer::calTrkIso(const reco::GsfElectron& ele,
0281 IsolationCalculators const& isolationCalculators) {
0282 if (ele.gsfTrack().isNull()) {
0283 return std::numeric_limits<float>::max();
0284 }
0285
0286 float trkIso = 0.;
0287 for (auto& calculator : isolationCalculators) {
0288 trkIso += (*calculator)(*ele.gsfTrack()).ptSum;
0289 }
0290 return trkIso;
0291 }
0292
0293 template <typename T>
0294 void ElectronHEEPIDValueMapProducer::writeValueMap(edm::Event& iEvent,
0295 const edm::Handle<edm::View<reco::GsfElectron>>& handle,
0296 const std::vector<T>& values,
0297 const std::string& label) {
0298 std::unique_ptr<edm::ValueMap<T>> valMap(new edm::ValueMap<T>());
0299 typename edm::ValueMap<T>::Filler filler(*valMap);
0300 filler.insert(handle, values.begin(), values.end());
0301 filler.fill();
0302 iEvent.put(std::move(valMap), label);
0303 }
0304
0305 void ElectronHEEPIDValueMapProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0306 edm::ParameterSetDescription desc;
0307 desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0308 desc.add<edm::InputTag>("ebRecHitsAOD", edm::InputTag("reducedEcalRecHitsEB"));
0309 desc.add<edm::InputTag>("eeRecHitsAOD", edm::InputTag("reducedEcalRecHitsEE"));
0310 desc.add<std::vector<edm::InputTag>>("candsAOD", {edm::InputTag("packedCandidates")});
0311 desc.add<std::vector<std::string>>("candVetosAOD", {"none"});
0312 desc.add<edm::InputTag>("elesAOD", edm::InputTag("gedGsfElectrons"));
0313
0314 desc.add<edm::InputTag>("ebRecHitsMiniAOD", edm::InputTag("reducedEcalRecHitsEB"));
0315 desc.add<edm::InputTag>("eeRecHitsMiniAOD", edm::InputTag("reducedEcalRecHitsEE"));
0316 desc.add<std::vector<edm::InputTag>>("candsMiniAOD", {edm::InputTag("packedCandidates")});
0317 desc.add<std::vector<std::string>>("candVetosMiniAOD", {"none"});
0318 desc.add<edm::InputTag>("elesMiniAOD", edm::InputTag("gedGsfElectrons"));
0319 desc.add<int>("dataFormat", 0);
0320 desc.add<bool>("makeTrkIso04", false);
0321 desc.add("trkIsoConfig", EleTkIsolFromCands::pSetDescript());
0322 desc.add("trkIso04Config", EleTkIsolFromCands::pSetDescript());
0323
0324 descriptions.addDefault(desc);
0325 }
0326
0327 DEFINE_FWK_MODULE(ElectronHEEPIDValueMapProducer);
