Skims/src/IsoPhotonEBSelector.cc

0001 // -*- C++ -*-
0002 // Class:      IsoPhotonEBSelector
0003 //
0004 // Original Author:  Silvia Taroni
0005 //         Created:  Wed, 29 Nov 2017 18:23:54 GMT
0006 //
0007 //
0008
0009 #include "FWCore/PluginManager/interface/ModuleDef.h"
0010
0011 // system include files
0012
0013 #include <algorithm>
0014 #include <iostream>
0015 #include <memory>
0016 #include <string>
0017 #include <vector>
0018
0019 // user include files
0020 #include "FWCore/Framework/interface/Frameworkfwd.h"
0021 #include "FWCore/Framework/interface/global/EDFilter.h"
0022
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0028
0029 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0030 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0031 //
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 #include "FWCore/Framework/interface/ConsumesCollector.h"
0034 #include "CommonTools/UtilAlgos/interface/SingleObjectSelector.h"
0035 #include "DataFormats/Common/interface/View.h"
0036
0037 using namespace std;
0038 using namespace reco;
0039 namespace edm {
0040   class EventSetup;
0041 }
0042
0043 class IsoPhotonEBSelector {
0044 public:
0045   IsoPhotonEBSelector(const edm::ParameterSet&, edm::ConsumesCollector& iC);
0046
0047   static void fillPSetDescription(edm::ParameterSetDescription& desc);
0048
0049   bool operator()(const reco::Photon&) const;
0050   void newEvent(const edm::Event&, const edm::EventSetup&);
0051   const float getEffectiveArea(float eta) const;
0052   void printEffectiveAreas() const;
0053
0054   const edm::EDGetTokenT<double> theRhoToken_;
0055   edm::Handle<double> rhoHandle_;
0056
0057   const std::vector<double> absEtaMin_;            // low limit of the eta range
0058   const std::vector<double> absEtaMax_;            // upper limit of the eta range
0059   const std::vector<double> effectiveAreaValues_;  // effective area for this eta range
0060
0061   const edm::ParameterSet phIDWP_;
0062
0063   const vector<double> sigmaIEtaIEtaCut_;
0064   const vector<double> hOverECut_;
0065   const vector<double> relCombIso_;
0066 };
0067
0068 void IsoPhotonEBSelector::printEffectiveAreas() const {
0069   printf("  eta_min   eta_max    effective area\n");
0070   uint nEtaBins = absEtaMin_.size();
0071   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0072     printf("  %8.4f    %8.4f   %8.5f\n", absEtaMin_[iEta], absEtaMax_[iEta], effectiveAreaValues_[iEta]);
0073   }
0074 }
0075 const float IsoPhotonEBSelector::getEffectiveArea(float eta) const {
0076   float effArea = 0;
0077   uint nEtaBins = absEtaMin_.size();
0078   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0079     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0080       effArea = effectiveAreaValues_[iEta];
0081       break;
0082     }
0083   }
0084
0085   return effArea;
0086 }
0087
0088 void IsoPhotonEBSelector::fillPSetDescription(edm::ParameterSetDescription& desc) {
0089   desc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetCentralCalo"));
0090   desc.add<std::vector<double>>("absEtaMin", {0.0000, 1.0000, 1.4790, 2.0000, 2.2000, 2.3000, 2.4000});
0091   desc.add<std::vector<double>>("absEtaMax", {1.0000, 1.4790, 2.0000, 2.2000, 2.3000, 2.4000, 5.0000});
0092   desc.add<std::vector<double>>("effectiveAreaValues", {0.1703, 0.1715, 0.1213, 0.1230, 0.1635, 0.1937, 0.2393});
0093
0094   // Define the photonIDWP PSet
0095   edm::ParameterSetDescription photonIDWP;
0096   photonIDWP.add<std::vector<double>>("full5x5_sigmaIEtaIEtaCut", {0.011, -1.0});
0097   photonIDWP.add<std::vector<double>>("hOverECut", {0.1, -1.0});
0098   photonIDWP.add<std::vector<double>>("relCombIsolationWithEACut", {0.05, -1.0});
0099
0100   // Add the PSet to the main description
0101   desc.add<edm::ParameterSetDescription>("phID", photonIDWP);
0102 }
0103
0104 IsoPhotonEBSelector::IsoPhotonEBSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0105     : theRhoToken_(iC.consumes<double>(cfg.getParameter<edm::InputTag>("rho"))),
0106       absEtaMin_(cfg.getParameter<std::vector<double>>("absEtaMin")),
0107       absEtaMax_(cfg.getParameter<std::vector<double>>("absEtaMax")),
0108       effectiveAreaValues_(cfg.getParameter<std::vector<double>>("effectiveAreaValues")),
0109       phIDWP_(cfg.getParameter<edm::ParameterSet>("phID")),
0110       sigmaIEtaIEtaCut_(phIDWP_.getParameter<std::vector<double>>("full5x5_sigmaIEtaIEtaCut")),
0111       hOverECut_(phIDWP_.getParameter<std::vector<double>>("hOverECut")),
0112       relCombIso_(phIDWP_.getParameter<std::vector<double>>("relCombIsolationWithEACut")) {
0113   //printEffectiveAreas();
0114 }
0115
0116 void IsoPhotonEBSelector::newEvent(const edm::Event& ev, const edm::EventSetup&) {
0117   ev.getByToken(theRhoToken_, rhoHandle_);
0118 }
0119
0120 bool IsoPhotonEBSelector::operator()(const reco::Photon& ph) const {
0121   float pt_e = ph.pt();
0122   unsigned int ind = 0;
0123   float abseta = fabs((ph.superCluster().get())->position().eta());
0124
0125   if (ph.isEB()) {
0126     if (abseta > 1.479)
0127       return false;  // check if it is really needed
0128   }
0129   if (ph.isEE()) {
0130     ind = 1;
0131     if (abseta <= 1.479 || abseta >= 2.5)
0132       return false;  // check if it is really needed
0133   }
0134
0135   if (ph.full5x5_sigmaIetaIeta() > sigmaIEtaIEtaCut_[ind])
0136     return false;
0137   if (ph.hadronicOverEm() > hOverECut_[ind])
0138     return false;
0139   const float eA = getEffectiveArea(abseta);
0140   const float rho = rhoHandle_.isValid() ? (float)(*rhoHandle_.product()) : 0;
0141   if ((ph.getPflowIsolationVariables().chargedHadronIso +
0142        std::max(float(0.0),
0143                 ph.getPflowIsolationVariables().neutralHadronIso + ph.getPflowIsolationVariables().photonIso -
0144                     eA * rho)) > relCombIso_[ind] * pt_e)
0145     return false;
0146
0147   return true;
0148 }
0149
0150 EVENTSETUP_STD_INIT(IsoPhotonEBSelector);
0151
0152 typedef SingleObjectSelector<edm::View<reco::Photon>, IsoPhotonEBSelector> IsoPhotonEBSelectorAndSkim;
0153
0154 DEFINE_FWK_MODULE(IsoPhotonEBSelectorAndSkim);