Skims/src/IsoPhotonEBSelector.cc

0001 // -*- C++ -*-
0002 // Class:      ZElectronsSelector
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
0028 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0029 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0030 //
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 #include "FWCore/Framework/interface/ConsumesCollector.h"
0033 #include "CommonTools/UtilAlgos/interface/SingleObjectSelector.h"
0034 #include "DataFormats/Common/interface/View.h"
0035
0036 using namespace std;
0037 using namespace reco;
0038 namespace edm {
0039   class EventSetup;
0040 }
0041
0042 class IsoPhotonEBSelector {
0043 public:
0044   IsoPhotonEBSelector(const edm::ParameterSet&, edm::ConsumesCollector& iC);
0045   bool operator()(const reco::Photon&) const;
0046   void newEvent(const edm::Event&, const edm::EventSetup&);
0047   const float getEffectiveArea(float eta) const;
0048   void printEffectiveAreas() const;
0049
0050   const edm::EDGetTokenT<double> theRhoToken_;
0051   edm::Handle<double> rhoHandle_;
0052
0053   const std::vector<double> absEtaMin_;            // low limit of the eta range
0054   const std::vector<double> absEtaMax_;            // upper limit of the eta range
0055   const std::vector<double> effectiveAreaValues_;  // effective area for this eta range
0056
0057   const edm::ParameterSet phIDWP_;
0058
0059   const vector<double> sigmaIEtaIEtaCut_;
0060   const vector<double> hOverECut_;
0061   const vector<double> relCombIso_;
0062 };
0063
0064 void IsoPhotonEBSelector::printEffectiveAreas() const {
0065   printf("  eta_min   eta_max    effective area\n");
0066   uint nEtaBins = absEtaMin_.size();
0067   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0068     printf("  %8.4f    %8.4f   %8.5f\n", absEtaMin_[iEta], absEtaMax_[iEta], effectiveAreaValues_[iEta]);
0069   }
0070 }
0071 const float IsoPhotonEBSelector::getEffectiveArea(float eta) const {
0072   float effArea = 0;
0073   uint nEtaBins = absEtaMin_.size();
0074   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0075     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0076       effArea = effectiveAreaValues_[iEta];
0077       break;
0078     }
0079   }
0080
0081   return effArea;
0082 }
0083
0084 IsoPhotonEBSelector::IsoPhotonEBSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0085     : theRhoToken_(iC.consumes<double>(cfg.getParameter<edm::InputTag>("rho"))),
0086       absEtaMin_(cfg.getParameter<std::vector<double> >("absEtaMin")),
0087       absEtaMax_(cfg.getParameter<std::vector<double> >("absEtaMax")),
0088       effectiveAreaValues_(cfg.getParameter<std::vector<double> >("effectiveAreaValues")),
0089       phIDWP_(cfg.getParameter<edm::ParameterSet>("phID")),
0090       sigmaIEtaIEtaCut_(phIDWP_.getParameter<std::vector<double> >("full5x5_sigmaIEtaIEtaCut")),
0091       hOverECut_(phIDWP_.getParameter<std::vector<double> >("hOverECut")),
0092       relCombIso_(phIDWP_.getParameter<std::vector<double> >("relCombIsolationWithEACut")) {
0093   //printEffectiveAreas();
0094 }
0095
0096 void IsoPhotonEBSelector::newEvent(const edm::Event& ev, const edm::EventSetup&) {
0097   ev.getByToken(theRhoToken_, rhoHandle_);
0098 }
0099
0100 bool IsoPhotonEBSelector::operator()(const reco::Photon& ph) const {
0101   float pt_e = ph.pt();
0102   unsigned int ind = 0;
0103   float abseta = fabs((ph.superCluster().get())->position().eta());
0104
0105   if (ph.isEB()) {
0106     if (abseta > 1.479)
0107       return false;  // check if it is really needed
0108   }
0109   if (ph.isEE()) {
0110     ind = 1;
0111     if (abseta <= 1.479 || abseta >= 2.5)
0112       return false;  // check if it is really needed
0113   }
0114
0115   if (ph.full5x5_sigmaIetaIeta() > sigmaIEtaIEtaCut_[ind])
0116     return false;
0117   if (ph.hadronicOverEm() > hOverECut_[ind])
0118     return false;
0119   const float eA = getEffectiveArea(abseta);
0120   const float rho = rhoHandle_.isValid() ? (float)(*rhoHandle_.product()) : 0;
0121   if ((ph.getPflowIsolationVariables().chargedHadronIso +
0122        std::max(float(0.0),
0123                 ph.getPflowIsolationVariables().neutralHadronIso + ph.getPflowIsolationVariables().photonIso -
0124                     eA * rho)) > relCombIso_[ind] * pt_e)
0125     return false;
0126
0127   return true;
0128 }
0129
0130 EVENTSETUP_STD_INIT(IsoPhotonEBSelector);
0131
0132 typedef SingleObjectSelector<edm::View<reco::Photon>, IsoPhotonEBSelector> IsoPhotonEBSelectorAndSkim;
0133
0134 DEFINE_FWK_MODULE(IsoPhotonEBSelectorAndSkim);