Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:41 

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         doQuadratic_(iConfig.getParameter<bool>("doQuadratic")) {
0050     if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0051         typeid(T) == typeid(pat::IsolatedTrack)) {
0052       produces<edm::ValueMap<float>>("miniIsoChg");
0053       produces<edm::ValueMap<float>>("miniIsoAll");
0054       ea_miniiso_ =
0055           std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_MiniIso")).fullPath());
0056       rho_miniiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_MiniIso"));
0057     }
0058     if ((typeid(T) == typeid(pat::Electron))) {
0059       produces<edm::ValueMap<float>>("PFIsoChg");
0060       produces<edm::ValueMap<float>>("PFIsoAll");
0061       produces<edm::ValueMap<float>>("PFIsoAll04");
0062       ea_pfiso_ = std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso")).fullPath());
0063       rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
0064     } else if ((typeid(T) == typeid(pat::Photon))) {
0065       rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
0066 
0067       if (!doQuadratic_) {
0068         produces<edm::ValueMap<float>>("PFIsoChg");
0069         produces<edm::ValueMap<float>>("PFIsoAll");
0070 
0071         ea_pfiso_chg_ =
0072             std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Chg")).fullPath());
0073         ea_pfiso_neu_ =
0074             std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Neu")).fullPath());
0075         ea_pfiso_pho_ =
0076             std::make_unique<EffectiveAreas>((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Pho")).fullPath());
0077 
0078       }
0079 
0080       else {
0081         produces<edm::ValueMap<float>>("PFIsoChgQuadratic");
0082         produces<edm::ValueMap<float>>("PFIsoAllQuadratic");
0083 
0084         quadratic_ea_pfiso_chg_ = std::make_unique<EffectiveAreas>(
0085             (iConfig.getParameter<edm::FileInPath>("QuadraticEAFile_PFIso_Chg")).fullPath(), true);
0086         quadratic_ea_pfiso_ecal_ = std::make_unique<EffectiveAreas>(
0087             (iConfig.getParameter<edm::FileInPath>("QuadraticEAFile_PFIso_ECal")).fullPath(), true);
0088         quadratic_ea_pfiso_hcal_ = std::make_unique<EffectiveAreas>(
0089             (iConfig.getParameter<edm::FileInPath>("QuadraticEAFile_PFIso_HCal")).fullPath(), true);
0090       }
0091     }
0092   }
0093 
0094   ~IsoValueMapProducer() override {}
0095 
0096   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0097 
0098 private:
0099   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0100 
0101   // ----------member data ---------------------------
0102 
0103   edm::EDGetTokenT<edm::View<T>> src_;
0104   bool relative_;
0105   bool doQuadratic_;
0106   edm::EDGetTokenT<double> rho_miniiso_;
0107   edm::EDGetTokenT<double> rho_pfiso_;
0108   std::unique_ptr<EffectiveAreas> ea_miniiso_;
0109   std::unique_ptr<EffectiveAreas> ea_pfiso_;
0110   std::unique_ptr<EffectiveAreas> ea_pfiso_chg_;
0111   std::unique_ptr<EffectiveAreas> ea_pfiso_neu_;
0112   std::unique_ptr<EffectiveAreas> ea_pfiso_pho_;
0113   std::unique_ptr<EffectiveAreas> quadratic_ea_pfiso_chg_;
0114   std::unique_ptr<EffectiveAreas> quadratic_ea_pfiso_ecal_;
0115   std::unique_ptr<EffectiveAreas> quadratic_ea_pfiso_hcal_;
0116 
0117   float getEtaForEA(const T*) const;
0118   void doMiniIso(edm::Event&) const;
0119   void doPFIsoEle(edm::Event&) const;
0120   void doPFIsoPho(edm::Event&) const;
0121   void doPFIsoPhoQuadratic(edm::Event&) const;
0122 };
0123 
0124 //
0125 // constants, enums and typedefs
0126 //
0127 
0128 //
0129 // static data member definitions
0130 //
0131 
0132 template <typename T>
0133 float IsoValueMapProducer<T>::getEtaForEA(const T* obj) const {
0134   return obj->eta();
0135 }
0136 template <>
0137 float IsoValueMapProducer<pat::Electron>::getEtaForEA(const pat::Electron* el) const {
0138   return el->superCluster()->eta();
0139 }
0140 template <>
0141 float IsoValueMapProducer<pat::Photon>::getEtaForEA(const pat::Photon* ph) const {
0142   return ph->superCluster()->eta();
0143 }
0144 
0145 template <typename T>
0146 void IsoValueMapProducer<T>::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0147   if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0148       typeid(T) == typeid(pat::IsolatedTrack)) {
0149     doMiniIso(iEvent);
0150   };
0151   if ((typeid(T) == typeid(pat::Electron))) {
0152     doPFIsoEle(iEvent);
0153   }
0154   if ((typeid(T) == typeid(pat::Photon))) {
0155     if (!doQuadratic_)
0156       doPFIsoPho(iEvent);
0157     else
0158       doPFIsoPhoQuadratic(iEvent);
0159   }
0160 }
0161 
0162 template <typename T>
0163 void IsoValueMapProducer<T>::doMiniIso(edm::Event& iEvent) const {
0164   auto src = iEvent.getHandle(src_);
0165   const auto& rho = iEvent.get(rho_miniiso_);
0166 
0167   unsigned int nInput = src->size();
0168 
0169   std::vector<float> miniIsoChg, miniIsoAll;
0170   miniIsoChg.reserve(nInput);
0171   miniIsoAll.reserve(nInput);
0172 
0173   for (const auto& obj : *src) {
0174     auto iso = obj.miniPFIsolation();
0175     auto chg = iso.chargedHadronIso();
0176     auto neu = iso.neutralHadronIso();
0177     auto pho = iso.photonIso();
0178     auto ea = ea_miniiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0179     float R = 10.0 / std::min(std::max(obj.pt(), 50.0), 200.0);
0180     ea *= std::pow(R / 0.3, 2);
0181     float scale = relative_ ? 1.0 / obj.pt() : 1;
0182     miniIsoChg.push_back(scale * chg);
0183     miniIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - rho * ea)));
0184   }
0185 
0186   auto miniIsoChgV = std::make_unique<edm::ValueMap<float>>();
0187   edm::ValueMap<float>::Filler fillerChg(*miniIsoChgV);
0188   fillerChg.insert(src, miniIsoChg.begin(), miniIsoChg.end());
0189   fillerChg.fill();
0190   auto miniIsoAllV = std::make_unique<edm::ValueMap<float>>();
0191   edm::ValueMap<float>::Filler fillerAll(*miniIsoAllV);
0192   fillerAll.insert(src, miniIsoAll.begin(), miniIsoAll.end());
0193   fillerAll.fill();
0194 
0195   iEvent.put(std::move(miniIsoChgV), "miniIsoChg");
0196   iEvent.put(std::move(miniIsoAllV), "miniIsoAll");
0197 }
0198 
0199 template <>
0200 void IsoValueMapProducer<pat::Photon>::doMiniIso(edm::Event& iEvent) const {}
0201 
0202 template <typename T>
0203 void IsoValueMapProducer<T>::doPFIsoEle(edm::Event& iEvent) const {}
0204 
0205 template <>
0206 void IsoValueMapProducer<pat::Electron>::doPFIsoEle(edm::Event& iEvent) const {
0207   edm::Handle<edm::View<pat::Electron>> src;
0208   iEvent.getByToken(src_, src);
0209   const auto& rho = iEvent.get(rho_pfiso_);
0210 
0211   unsigned int nInput = src->size();
0212 
0213   std::vector<float> PFIsoChg, PFIsoAll, PFIsoAll04;
0214   PFIsoChg.reserve(nInput);
0215   PFIsoAll.reserve(nInput);
0216   PFIsoAll04.reserve(nInput);
0217 
0218   for (const auto& obj : *src) {
0219     auto iso = obj.pfIsolationVariables();
0220     auto chg = iso.sumChargedHadronPt;
0221     auto neu = iso.sumNeutralHadronEt;
0222     auto pho = iso.sumPhotonEt;
0223     auto ea = ea_pfiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0224     float scale = relative_ ? 1.0 / obj.pt() : 1;
0225     PFIsoChg.push_back(scale * chg);
0226     PFIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - rho * ea)));
0227     PFIsoAll04.push_back(scale * (obj.chargedHadronIso() +
0228                                   std::max(0.0, obj.neutralHadronIso() + obj.photonIso() - rho * ea * 16. / 9.)));
0229   }
0230 
0231   auto PFIsoChgV = std::make_unique<edm::ValueMap<float>>();
0232   edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
0233   fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
0234   fillerChg.fill();
0235   auto PFIsoAllV = std::make_unique<edm::ValueMap<float>>();
0236   edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
0237   fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
0238   fillerAll.fill();
0239   auto PFIsoAll04V = std::make_unique<edm::ValueMap<float>>();
0240   edm::ValueMap<float>::Filler fillerAll04(*PFIsoAll04V);
0241   fillerAll04.insert(src, PFIsoAll04.begin(), PFIsoAll04.end());
0242   fillerAll04.fill();
0243 
0244   iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
0245   iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
0246   iEvent.put(std::move(PFIsoAll04V), "PFIsoAll04");
0247 }
0248 
0249 template <typename T>
0250 void IsoValueMapProducer<T>::doPFIsoPho(edm::Event& iEvent) const {}
0251 
0252 template <>
0253 void IsoValueMapProducer<pat::Photon>::doPFIsoPho(edm::Event& iEvent) const {
0254   edm::Handle<edm::View<pat::Photon>> src;
0255   iEvent.getByToken(src_, src);
0256   const auto& rho = iEvent.get(rho_pfiso_);
0257 
0258   unsigned int nInput = src->size();
0259 
0260   std::vector<float> PFIsoChg, PFIsoAll;
0261 
0262   PFIsoChg.reserve(nInput);
0263   PFIsoAll.reserve(nInput);
0264 
0265   for (const auto& obj : *src) {
0266     auto chg = obj.chargedHadronIso();
0267     auto neu = obj.neutralHadronIso();
0268     auto pho = obj.photonIso();
0269 
0270     auto ea_chg = ea_pfiso_chg_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0271     auto ea_neu = ea_pfiso_neu_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0272     auto ea_pho = ea_pfiso_pho_->getEffectiveArea(fabs(getEtaForEA(&obj)));
0273 
0274     float scale = relative_ ? 1.0 / obj.pt() : 1;
0275     PFIsoChg.push_back(scale * std::max(0.0, chg - rho * ea_chg));
0276     PFIsoAll.push_back(PFIsoChg.back() +
0277                        scale * (std::max(0.0, neu - rho * ea_neu) + std::max(0.0, pho - rho * ea_pho)));
0278   }
0279 
0280   auto PFIsoChgV = std::make_unique<edm::ValueMap<float>>();
0281   edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
0282   fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
0283   fillerChg.fill();
0284   auto PFIsoAllV = std::make_unique<edm::ValueMap<float>>();
0285   edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
0286   fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
0287   fillerAll.fill();
0288 
0289   iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
0290   iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
0291 }
0292 
0293 template <typename T>
0294 void IsoValueMapProducer<T>::doPFIsoPhoQuadratic(edm::Event& iEvent) const {}
0295 
0296 template <>
0297 void IsoValueMapProducer<pat::Photon>::doPFIsoPhoQuadratic(edm::Event& iEvent) const {
0298   edm::Handle<edm::View<pat::Photon>> src;
0299   iEvent.getByToken(src_, src);
0300   const auto& rho = iEvent.get(rho_pfiso_);
0301 
0302   unsigned int nInput = src->size();
0303 
0304   std::vector<float> PFIsoChgQuadratic, PFIsoAllQuadratic;
0305 
0306   PFIsoChgQuadratic.reserve(nInput);
0307   PFIsoAllQuadratic.reserve(nInput);
0308 
0309   for (const auto& obj : *src) {
0310     auto chg = obj.chargedHadronIso();
0311     auto ecal = obj.ecalPFClusterIso();
0312     auto hcal = obj.hcalPFClusterIso();
0313 
0314     auto quadratic_ea_chg = quadratic_ea_pfiso_chg_->getQuadraticEA(fabs(getEtaForEA(&obj)));
0315     auto linear_ea_chg = quadratic_ea_pfiso_chg_->getLinearEA(fabs(getEtaForEA(&obj)));
0316     auto quadratic_ea_ecal = quadratic_ea_pfiso_ecal_->getQuadraticEA(fabs(getEtaForEA(&obj)));
0317     auto linear_ea_ecal = quadratic_ea_pfiso_ecal_->getLinearEA(fabs(getEtaForEA(&obj)));
0318     auto quadratic_ea_hcal = quadratic_ea_pfiso_hcal_->getQuadraticEA(fabs(getEtaForEA(&obj)));
0319     auto linear_ea_hcal = quadratic_ea_pfiso_hcal_->getLinearEA(fabs(getEtaForEA(&obj)));
0320 
0321     float scale = relative_ ? 1.0 / obj.pt() : 1;
0322 
0323     PFIsoChgQuadratic.push_back(scale * std::max(0.0, chg - (quadratic_ea_chg * rho * rho + linear_ea_chg * rho)));
0324     PFIsoAllQuadratic.push_back(PFIsoChgQuadratic.back() +
0325                                 scale * (std::max(0.0, ecal - (quadratic_ea_ecal * rho * rho + linear_ea_ecal * rho)) +
0326                                          std::max(0.0, hcal - (quadratic_ea_hcal * rho * rho + linear_ea_hcal * rho))));
0327   }
0328 
0329   auto PFIsoChgQuadraticV = std::make_unique<edm::ValueMap<float>>();
0330   edm::ValueMap<float>::Filler fillerChgQuadratic(*PFIsoChgQuadraticV);
0331   fillerChgQuadratic.insert(src, PFIsoChgQuadratic.begin(), PFIsoChgQuadratic.end());
0332   fillerChgQuadratic.fill();
0333 
0334   auto PFIsoAllQuadraticV = std::make_unique<edm::ValueMap<float>>();
0335   edm::ValueMap<float>::Filler fillerAllQuadratic(*PFIsoAllQuadraticV);
0336   fillerAllQuadratic.insert(src, PFIsoAllQuadratic.begin(), PFIsoAllQuadratic.end());
0337   fillerAllQuadratic.fill();
0338 
0339   iEvent.put(std::move(PFIsoChgQuadraticV), "PFIsoChgQuadratic");
0340   iEvent.put(std::move(PFIsoAllQuadraticV), "PFIsoAllQuadratic");
0341 }
0342 
0343 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0344 template <typename T>
0345 void IsoValueMapProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0346   edm::ParameterSetDescription desc;
0347   desc.add<edm::InputTag>("src")->setComment("input physics object collection");
0348   desc.add<bool>("relative")->setComment("compute relative isolation instead of absolute one");
0349   desc.add<bool>("doQuadratic", false)->setComment("flag to do quadratic EA corections for photons");
0350   if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
0351       typeid(T) == typeid(pat::IsolatedTrack)) {
0352     desc.add<edm::FileInPath>("EAFile_MiniIso")
0353         ->setComment("txt file containing effective areas to be used for mini-isolation pileup subtraction");
0354     desc.add<edm::InputTag>("rho_MiniIso")
0355         ->setComment("rho to be used for effective-area based mini-isolation pileup subtraction");
0356   }
0357   if ((typeid(T) == typeid(pat::Electron))) {
0358     desc.add<edm::FileInPath>("EAFile_PFIso")
0359         ->setComment(
0360             "txt file containing effective areas to be used for PF-isolation pileup subtraction for electrons");
0361     desc.add<edm::InputTag>("rho_PFIso")
0362         ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for electrons");
0363   }
0364   if ((typeid(T) == typeid(pat::Photon))) {
0365     desc.addOptional<edm::InputTag>("mapIsoChg")
0366         ->setComment("input charged PF isolation calculated in VID for photons");
0367     desc.addOptional<edm::InputTag>("mapIsoNeu")
0368         ->setComment("input neutral PF isolation calculated in VID for photons");
0369     desc.addOptional<edm::InputTag>("mapIsoPho")->setComment("input photon PF isolation calculated in VID for photons");
0370 
0371     desc.addOptional<edm::FileInPath>("EAFile_PFIso_Chg")
0372         ->setComment(
0373             "txt file containing effective areas to be used for charged PF-isolation pileup subtraction for photons");
0374     desc.addOptional<edm::FileInPath>("EAFile_PFIso_Neu")
0375         ->setComment(
0376             "txt file containing effective areas to be used for neutral PF-isolation pileup subtraction for photons");
0377     desc.addOptional<edm::FileInPath>("EAFile_PFIso_Pho")
0378         ->setComment(
0379             "txt file containing effective areas to be used for photon PF-isolation pileup subtraction for photons");
0380 
0381     desc.add<edm::InputTag>("rho_PFIso")
0382         ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for photons");
0383 
0384     desc.addOptional<edm::FileInPath>("QuadraticEAFile_PFIso_Chg")
0385         ->setComment(
0386             "txt file containing quadratic effective areas to be used for charged PF-isolation pileup subtraction for "
0387             "photons");
0388     desc.addOptional<edm::FileInPath>("QuadraticEAFile_PFIso_ECal")
0389         ->setComment(
0390             "txt file containing quadratic effective areas to be used for ecal PF-isolation pileup subtraction for "
0391             "photons");
0392     desc.addOptional<edm::FileInPath>("QuadraticEAFile_PFIso_HCal")
0393         ->setComment(
0394             "txt file containing quadratic effective areas to be used for hcal PF-isolation pileup subtraction for "
0395             "photons");
0396   }
0397 
0398   std::string modname;
0399   if (typeid(T) == typeid(pat::Muon))
0400     modname += "Muon";
0401   else if (typeid(T) == typeid(pat::Electron))
0402     modname += "Ele";
0403   else if (typeid(T) == typeid(pat::Photon))
0404     modname += "Pho";
0405   else if (typeid(T) == typeid(pat::IsolatedTrack))
0406     modname += "IsoTrack";
0407   modname += "IsoValueMapProducer";
0408   descriptions.add(modname, desc);
0409 }
0410 
0411 typedef IsoValueMapProducer<pat::Muon> MuonIsoValueMapProducer;
0412 typedef IsoValueMapProducer<pat::Electron> EleIsoValueMapProducer;
0413 typedef IsoValueMapProducer<pat::Photon> PhoIsoValueMapProducer;
0414 typedef IsoValueMapProducer<pat::IsolatedTrack> IsoTrackIsoValueMapProducer;
0415 
0416 //define this as a plug-in
0417 DEFINE_FWK_MODULE(MuonIsoValueMapProducer);
0418 DEFINE_FWK_MODULE(EleIsoValueMapProducer);
0419 DEFINE_FWK_MODULE(PhoIsoValueMapProducer);
0420 DEFINE_FWK_MODULE(IsoTrackIsoValueMapProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!