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