Skims/src/ZElectronsSelector.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/GsfElectron.h"
0029 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0030 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0031 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0032 #include "CommonTools/UtilAlgos/interface/SingleObjectSelector.h"
0033 //
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035 #include "FWCore/Framework/interface/ConsumesCollector.h"
0036 #include "CommonTools/UtilAlgos/interface/SingleObjectSelector.h"
0037 #include "DataFormats/Common/interface/View.h"
0038
0039 using namespace std;
0040 using namespace reco;
0041 namespace edm {
0042   class EventSetup;
0043 }
0044
0045 class ZElectronsSelector {
0046 public:
0047   ZElectronsSelector(const edm::ParameterSet&, edm::ConsumesCollector& iC);
0048   bool operator()(const reco::GsfElectron&) const;
0049   void newEvent(const edm::Event&, const edm::EventSetup&);
0050   const float getEffectiveArea(float eta) const;
0051   void printEffectiveAreas() const;
0052
0053   edm::EDGetTokenT<double> theRhoToken;
0054   edm::EDGetTokenT<reco::GsfElectronCollection> theGsfEToken;
0055   edm::Handle<double> _rhoHandle;
0056
0057   std::vector<double> absEtaMin_;            // low limit of the eta range
0058   std::vector<double> absEtaMax_;            // upper limit of the eta range
0059   std::vector<double> effectiveAreaValues_;  // effective area for this eta range
0060
0061   edm::ParameterSet eleIDWP;
0062
0063   vector<int> missHits;
0064   vector<double> sigmaIEtaIEtaCut;
0065   vector<double> dEtaInSeedCut;
0066   vector<double> dPhiInCut;
0067   vector<double> hOverECut;
0068   vector<double> relCombIso;
0069   vector<double> EInvMinusPInv;
0070 };
0071
0072 void ZElectronsSelector::printEffectiveAreas() const {
0073   printf("  eta_min   eta_max    effective area\n");
0074   uint nEtaBins = absEtaMin_.size();
0075   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0076     printf("  %8.4f    %8.4f   %8.5f\n", absEtaMin_[iEta], absEtaMax_[iEta], effectiveAreaValues_[iEta]);
0077   }
0078 }
0079 const float ZElectronsSelector::getEffectiveArea(float eta) const {
0080   float effArea = 0;
0081   uint nEtaBins = absEtaMin_.size();
0082   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0083     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0084       effArea = effectiveAreaValues_[iEta];
0085       break;
0086     }
0087   }
0088
0089   return effArea;
0090 }
0091
0092 ZElectronsSelector::ZElectronsSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0093     : theRhoToken(iC.consumes<double>(cfg.getParameter<edm::InputTag>("rho"))) {
0094   absEtaMin_ = cfg.getParameter<std::vector<double> >("absEtaMin");
0095   absEtaMax_ = cfg.getParameter<std::vector<double> >("absEtaMax");
0096   effectiveAreaValues_ = cfg.getParameter<std::vector<double> >("effectiveAreaValues");
0097   //printEffectiveAreas();
0098   eleIDWP = cfg.getParameter<edm::ParameterSet>("eleID");
0099
0100   missHits = eleIDWP.getParameter<std::vector<int> >("missingHitsCut");
0101   sigmaIEtaIEtaCut = eleIDWP.getParameter<std::vector<double> >("full5x5_sigmaIEtaIEtaCut");
0102   dEtaInSeedCut = eleIDWP.getParameter<std::vector<double> >("dEtaInSeedCut");
0103   dPhiInCut = eleIDWP.getParameter<std::vector<double> >("dPhiInCut");
0104   hOverECut = eleIDWP.getParameter<std::vector<double> >("hOverECut");
0105   relCombIso = eleIDWP.getParameter<std::vector<double> >("relCombIsolationWithEACut");
0106   EInvMinusPInv = eleIDWP.getParameter<std::vector<double> >("EInverseMinusPInverseCut");
0107 }
0108
0109 void ZElectronsSelector::newEvent(const edm::Event& ev, const edm::EventSetup&) {
0110   ev.getByToken(theRhoToken, _rhoHandle);
0111 }
0112
0113 bool ZElectronsSelector::operator()(const reco::GsfElectron& el) const {
0114   float pt_e = el.pt();
0115   unsigned int ind = 0;
0116   auto etrack = el.gsfTrack();
0117   float abseta = fabs((el.superCluster().get())->position().eta());
0118
0119   if (el.isEB()) {
0120     if (abseta > 1.479)
0121       return false;  // check if it is really needed
0122   }
0123   if (el.isEE()) {
0124     ind = 1;
0125     if (abseta < 1.479)
0126       return false;  // check if it is really needed
0127     if (abseta >= 2.5)
0128       return false;  // check if it is really needed
0129   }
0130
0131   if (etrack->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) > missHits[ind])
0132     return false;
0133   if (el.full5x5_sigmaIetaIeta() > sigmaIEtaIEtaCut[ind])
0134     return false;
0135   if (fabs(el.deltaPhiSuperClusterTrackAtVtx()) > dPhiInCut[ind])
0136     return false;
0137   if (fabs(el.deltaEtaSeedClusterTrackAtVtx()) > dEtaInSeedCut[ind])
0138     return false;
0139   if (el.hadronicOverEm() > hOverECut[ind])
0140     return false;
0141   const float eA = getEffectiveArea(abseta);
0142   const float rho = _rhoHandle.isValid() ? (float)(*_rhoHandle.product()) : 0;
0143   if ((el.pfIsolationVariables().sumChargedHadronPt +
0144        std::max(float(0.0),
0145                 el.pfIsolationVariables().sumNeutralHadronEt + el.pfIsolationVariables().sumPhotonEt - eA * rho)) >
0146       relCombIso[ind] * pt_e)
0147     return false;
0148   const float ecal_energy_inverse = 1.0 / el.ecalEnergy();
0149   const float eSCoverP = el.eSuperClusterOverP();
0150   if (std::abs(1.0 - eSCoverP) * ecal_energy_inverse > EInvMinusPInv[ind])
0151     return false;
0152
0153   return true;
0154 }
0155
0156 EVENTSETUP_STD_INIT(ZElectronsSelector);
0157
0158 typedef SingleObjectSelector<edm::View<reco::GsfElectron>, ZElectronsSelector> ZElectronsSelectorAndSkim;
0159
0160 DEFINE_FWK_MODULE(ZElectronsSelectorAndSkim);