Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-26 00:02:52

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 //Heavily inspired from ElectronIDValueMapProducer
0033 class ElectronHEEPIDValueMapProducer : public edm::stream::EDProducer<> {
0034 private:
0035   //helper classes to handle AOD vs MiniAOD
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);