Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:32:45

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/NanoAOD
0004 // Class:      IsoValueMapProducer
0005 //
0006 /**\class IsoValueMapProducer IsoValueMapProducer.cc PhysicsTools/NanoAOD/plugins/IsoValueMapProducer.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Marco Peruzzi
0015 //         Created:  Mon, 04 Sep 2017 22:43:53 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/global/EDProducer.h"
0025 
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031 
0032 #include "CommonTools/Egamma/interface/EffectiveAreas.h"
0033 
0034 #include "DataFormats/PatCandidates/interface/Muon.h"
0035 #include "DataFormats/PatCandidates/interface/Electron.h"
0036 #include "DataFormats/PatCandidates/interface/Photon.h"
0037 #include "DataFormats/PatCandidates/interface/IsolatedTrack.h"
0038 
0039 //
0040 // class declaration
0041 //
0042 
0043 template <typename T>
0044 class IsoValueMapProducer : public edm::global::EDProducer<> {
0045 public:
0046   explicit IsoValueMapProducer(const edm::ParameterSet& iConfig)
0047       : src_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("src"))),
0048         relative_(iConfig.getParameter<bool>("relative")) {
0049     if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0050         typeid(T) == typeid(pat::IsolatedTrack)) {
0051       produces<edm::ValueMap<float>>("miniIsoChg");
0052       produces<edm::ValueMap<float>>("miniIsoAll");
0053       ea_miniiso_ =
0054           std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_MiniIso")).fullPath());
0055       rho_miniiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_MiniIso"));
0056     }
0057     if ((typeid(T) == typeid(pat::Electron))) {
0058       produces<edm::ValueMap<float>>("PFIsoChg");
0059       produces<edm::ValueMap<float>>("PFIsoAll");
0060       produces<edm::ValueMap<float>>("PFIsoAll04");
0061       ea_pfiso_ = std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso")).fullPath());
0062       rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
0063     } else if ((typeid(T) == typeid(pat::Photon))) {
0064       produces<edm::ValueMap<float>>("PFIsoChg");
0065       produces<edm::ValueMap<float>>("PFIsoAll");
0066       ea_pfiso_chg_ =
0067           std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Chg")).fullPath());
0068       ea_pfiso_neu_ =
0069           std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Neu")).fullPath());
0070       ea_pfiso_pho_ =
0071           std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Pho")).fullPath());
0072       rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
0073     }
0074   }
0075   ~IsoValueMapProducer() override {}
0076 
0077   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0078 
0079 private:
0080   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0081 
0082   // ----------member data ---------------------------
0083 
0084   edm::EDGetTokenT<edm::View<T>> src_;
0085   bool relative_;
0086   edm::EDGetTokenT<double> rho_miniiso_;
0087   edm::EDGetTokenT<double> rho_pfiso_;
0088   std::unique_ptr<EffectiveAreas> ea_miniiso_;
0089   std::unique_ptr<EffectiveAreas> ea_pfiso_;
0090   std::unique_ptr<EffectiveAreas> ea_pfiso_chg_;
0091   std::unique_ptr<EffectiveAreas> ea_pfiso_neu_;
0092   std::unique_ptr<EffectiveAreas> ea_pfiso_pho_;
0093   float getEtaForEA(const T*) const;
0094   void doMiniIso(edm::Event&) const;
0095   void doPFIsoEle(edm::Event&) const;
0096   void doPFIsoPho(edm::Event&) const;
0097 };
0098 
0099 //
0100 // constants, enums and typedefs
0101 //
0102 
0103 //
0104 // static data member definitions
0105 //
0106 
0107 template <typename T>
0108 float IsoValueMapProducer<T>::getEtaForEA(const T* obj) const {
0109   return obj->eta();
0110 }
0111 template <>
0112 float IsoValueMapProducer<pat::Electron>::getEtaForEA(const pat::Electron* el) const {
0113   return el->superCluster()->eta();
0114 }
0115 template <>
0116 float IsoValueMapProducer<pat::Photon>::getEtaForEA(const pat::Photon* ph) const {
0117   return ph->superCluster()->eta();
0118 }
0119 
0120 template <typename T>
0121 void IsoValueMapProducer<T>::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0122   if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0123       typeid(T) == typeid(pat::IsolatedTrack)) {
0124     doMiniIso(iEvent);
0125   };
0126   if ((typeid(T) == typeid(pat::Electron))) {
0127     doPFIsoEle(iEvent);
0128   }
0129   if ((typeid(T) == typeid(pat::Photon))) {
0130     doPFIsoPho(iEvent);
0131   }
0132 }
0133 
0134 template <typename T>
0135 void IsoValueMapProducer<T>::doMiniIso(edm::Event& iEvent) const {
0136   edm::Handle<edm::View<T>> src;
0137   iEvent.getByToken(src_, src);
0138   edm::Handle<double> rho;
0139   iEvent.getByToken(rho_miniiso_, rho);
0140 
0141   unsigned int nInput = src->size();
0142 
0143   std::vector<float> miniIsoChg, miniIsoAll;
0144   miniIsoChg.reserve(nInput);
0145   miniIsoAll.reserve(nInput);
0146 
0147   for (const auto& obj : *src) {
0148     auto iso = obj.miniPFIsolation();
0149     auto chg = iso.chargedHadronIso();
0150     auto neu = iso.neutralHadronIso();
0151     auto pho = iso.photonIso();
0152     auto ea = ea_miniiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0153     float R = 10.0 / std::min(std::max(obj.pt(), 50.0), 200.0);
0154     ea *= std::pow(R / 0.3, 2);
0155     float scale = relative_ ? 1.0 / obj.pt() : 1;
0156     miniIsoChg.push_back(scale * chg);
0157     miniIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - (*rho) * ea)));
0158   }
0159 
0160   std::unique_ptr<edm::ValueMap<float>> miniIsoChgV(new edm::ValueMap<float>());
0161   edm::ValueMap<float>::Filler fillerChg(*miniIsoChgV);
0162   fillerChg.insert(src, miniIsoChg.begin(), miniIsoChg.end());
0163   fillerChg.fill();
0164   std::unique_ptr<edm::ValueMap<float>> miniIsoAllV(new edm::ValueMap<float>());
0165   edm::ValueMap<float>::Filler fillerAll(*miniIsoAllV);
0166   fillerAll.insert(src, miniIsoAll.begin(), miniIsoAll.end());
0167   fillerAll.fill();
0168 
0169   iEvent.put(std::move(miniIsoChgV), "miniIsoChg");
0170   iEvent.put(std::move(miniIsoAllV), "miniIsoAll");
0171 }
0172 
0173 template <>
0174 void IsoValueMapProducer<pat::Photon>::doMiniIso(edm::Event& iEvent) const {}
0175 
0176 template <typename T>
0177 void IsoValueMapProducer<T>::doPFIsoEle(edm::Event& iEvent) const {}
0178 
0179 template <>
0180 void IsoValueMapProducer<pat::Electron>::doPFIsoEle(edm::Event& iEvent) const {
0181   edm::Handle<edm::View<pat::Electron>> src;
0182   iEvent.getByToken(src_, src);
0183   edm::Handle<double> rho;
0184   iEvent.getByToken(rho_pfiso_, rho);
0185 
0186   unsigned int nInput = src->size();
0187 
0188   std::vector<float> PFIsoChg, PFIsoAll, PFIsoAll04;
0189   PFIsoChg.reserve(nInput);
0190   PFIsoAll.reserve(nInput);
0191   PFIsoAll04.reserve(nInput);
0192 
0193   for (const auto& obj : *src) {
0194     auto iso = obj.pfIsolationVariables();
0195     auto chg = iso.sumChargedHadronPt;
0196     auto neu = iso.sumNeutralHadronEt;
0197     auto pho = iso.sumPhotonEt;
0198     auto ea = ea_pfiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0199     float scale = relative_ ? 1.0 / obj.pt() : 1;
0200     PFIsoChg.push_back(scale * chg);
0201     PFIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - (*rho) * ea)));
0202     PFIsoAll04.push_back(scale * (obj.chargedHadronIso() +
0203                                   std::max(0.0, obj.neutralHadronIso() + obj.photonIso() - (*rho) * ea * 16. / 9.)));
0204   }
0205 
0206   std::unique_ptr<edm::ValueMap<float>> PFIsoChgV(new edm::ValueMap<float>());
0207   edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
0208   fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
0209   fillerChg.fill();
0210   std::unique_ptr<edm::ValueMap<float>> PFIsoAllV(new edm::ValueMap<float>());
0211   edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
0212   fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
0213   fillerAll.fill();
0214   std::unique_ptr<edm::ValueMap<float>> PFIsoAll04V(new edm::ValueMap<float>());
0215   edm::ValueMap<float>::Filler fillerAll04(*PFIsoAll04V);
0216   fillerAll04.insert(src, PFIsoAll04.begin(), PFIsoAll04.end());
0217   fillerAll04.fill();
0218 
0219   iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
0220   iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
0221   iEvent.put(std::move(PFIsoAll04V), "PFIsoAll04");
0222 }
0223 
0224 template <typename T>
0225 void IsoValueMapProducer<T>::doPFIsoPho(edm::Event& iEvent) const {}
0226 
0227 template <>
0228 void IsoValueMapProducer<pat::Photon>::doPFIsoPho(edm::Event& iEvent) const {
0229   edm::Handle<edm::View<pat::Photon>> src;
0230   iEvent.getByToken(src_, src);
0231   edm::Handle<double> rho;
0232   iEvent.getByToken(rho_pfiso_, rho);
0233 
0234   unsigned int nInput = src->size();
0235 
0236   std::vector<float> PFIsoChg, PFIsoAll;
0237   PFIsoChg.reserve(nInput);
0238   PFIsoAll.reserve(nInput);
0239 
0240   for (const auto& obj : *src) {
0241     auto chg = obj.chargedHadronIso();
0242     auto neu = obj.neutralHadronIso();
0243     auto pho = obj.photonIso();
0244     auto ea_chg = ea_pfiso_chg_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0245     auto ea_neu = ea_pfiso_neu_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0246     auto ea_pho = ea_pfiso_pho_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0247     float scale = relative_ ? 1.0 / obj.pt() : 1;
0248     PFIsoChg.push_back(scale * std::max(0.0, chg - (*rho) * ea_chg));
0249     PFIsoAll.push_back(PFIsoChg.back() +
0250                        scale * (std::max(0.0, neu - (*rho) * ea_neu) + std::max(0.0, pho - (*rho) * ea_pho)));
0251   }
0252 
0253   std::unique_ptr<edm::ValueMap<float>> PFIsoChgV(new edm::ValueMap<float>());
0254   edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
0255   fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
0256   fillerChg.fill();
0257   std::unique_ptr<edm::ValueMap<float>> PFIsoAllV(new edm::ValueMap<float>());
0258   edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
0259   fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
0260   fillerAll.fill();
0261 
0262   iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
0263   iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
0264 }
0265 
0266 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0267 template <typename T>
0268 void IsoValueMapProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   edm::ParameterSetDescription desc;
0270   desc.add<edm::InputTag>("src")->setComment("input physics object collection");
0271   desc.add<bool>("relative")->setComment("compute relative isolation instead of absolute one");
0272   if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0273       typeid(T) == typeid(pat::IsolatedTrack)) {
0274     desc.add<edm::FileInPath>("EAFile_MiniIso")
0275         ->setComment("txt file containing effective areas to be used for mini-isolation pileup subtraction");
0276     desc.add<edm::InputTag>("rho_MiniIso")
0277         ->setComment("rho to be used for effective-area based mini-isolation pileup subtraction");
0278   }
0279   if ((typeid(T) == typeid(pat::Electron))) {
0280     desc.add<edm::FileInPath>("EAFile_PFIso")
0281         ->setComment(
0282             "txt file containing effective areas to be used for PF-isolation pileup subtraction for electrons");
0283     desc.add<edm::InputTag>("rho_PFIso")
0284         ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for electrons");
0285   }
0286   if ((typeid(T) == typeid(pat::Photon))) {
0287     desc.add<edm::InputTag>("mapIsoChg")->setComment("input charged PF isolation calculated in VID for photons");
0288     desc.add<edm::InputTag>("mapIsoNeu")->setComment("input neutral PF isolation calculated in VID for photons");
0289     desc.add<edm::InputTag>("mapIsoPho")->setComment("input photon PF isolation calculated in VID for photons");
0290     desc.add<edm::FileInPath>("EAFile_PFIso_Chg")
0291         ->setComment(
0292             "txt file containing effective areas to be used for charged PF-isolation pileup subtraction for photons");
0293     desc.add<edm::FileInPath>("EAFile_PFIso_Neu")
0294         ->setComment(
0295             "txt file containing effective areas to be used for neutral PF-isolation pileup subtraction for photons");
0296     desc.add<edm::FileInPath>("EAFile_PFIso_Pho")
0297         ->setComment(
0298             "txt file containing effective areas to be used for photon PF-isolation pileup subtraction for photons");
0299     desc.add<edm::InputTag>("rho_PFIso")
0300         ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for photons");
0301   }
0302   std::string modname;
0303   if (typeid(T) == typeid(pat::Muon))
0304     modname += "Muon";
0305   else if (typeid(T) == typeid(pat::Electron))
0306     modname += "Ele";
0307   else if (typeid(T) == typeid(pat::Photon))
0308     modname += "Pho";
0309   else if (typeid(T) == typeid(pat::IsolatedTrack))
0310     modname += "IsoTrack";
0311   modname += "IsoValueMapProducer";
0312   descriptions.add(modname, desc);
0313 }
0314 
0315 typedef IsoValueMapProducer<pat::Muon> MuonIsoValueMapProducer;
0316 typedef IsoValueMapProducer<pat::Electron> EleIsoValueMapProducer;
0317 typedef IsoValueMapProducer<pat::Photon> PhoIsoValueMapProducer;
0318 typedef IsoValueMapProducer<pat::IsolatedTrack> IsoTrackIsoValueMapProducer;
0319 
0320 //define this as a plug-in
0321 DEFINE_FWK_MODULE(MuonIsoValueMapProducer);
0322 DEFINE_FWK_MODULE(EleIsoValueMapProducer);
0323 DEFINE_FWK_MODULE(PhoIsoValueMapProducer);
0324 DEFINE_FWK_MODULE(IsoTrackIsoValueMapProducer);