Lumi/src/ElectronIdentifier.cc

0001 #include "DQMOffline/Lumi/interface/ElectronIdentifier.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Utilities/interface/RegexMatch.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/VertexReco/interface/Vertex.h"
0006 #include "DataFormats/MuonReco/interface/Muon.h"
0007 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0008
0009 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0010 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0011 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0012 #include "CommonTools/Egamma/interface/ConversionTools.h"
0013
0014 #include <TLorentzVector.h>
0015 #include <TMath.h>
0016 #include <algorithm>
0017
0018 ElectronIdentifier::ElectronIdentifier(const edm::ParameterSet& c)
0019     : _effectiveAreas((c.getParameter<edm::FileInPath>("effAreasConfigFile")).fullPath())
0020
0021 {
0022   rho_ = -1;
0023   ID_ = -1;
0024   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.0115;
0025   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.011;
0026   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.00998;
0027   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.00998;
0028
0029   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.037;
0030   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.0314;
0031   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.0298;
0032   cuts_[EleIDCutNames::SIGMAIETA][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.0292;
0033
0034   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.00749;
0035   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.00477;
0036   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.00311;
0037   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.00308;
0038
0039   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.00895;
0040   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.00868;
0041   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.00609;
0042   cuts_[EleIDCutNames::DETAINSEED][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.00605;
0043
0044   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.228;
0045   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.222;
0046   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.103;
0047   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.0816;
0048
0049   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.213;
0050   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.213;
0051   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.045;
0052   cuts_[EleIDCutNames::DPHIIN][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.0394;
0053
0054   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.356;
0055   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.298;
0056   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.253;
0057   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.0414;
0058
0059   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.211;
0060   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.101;
0061   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.0878;
0062   cuts_[EleIDCutNames::HOVERE][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.0641;
0063
0064   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.175;
0065   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.0994;
0066   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.0695;
0067   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.0588;
0068
0069   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.159;
0070   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.107;
0071   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.0821;
0072   cuts_[EleIDCutNames::ISO][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.0571;
0073
0074   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 0.299;
0075   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 0.241;
0076   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 0.134;
0077   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 0.0129;
0078
0079   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 0.15;
0080   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 0.14;
0081   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 0.13;
0082   cuts_[EleIDCutNames::ONEOVERE][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 0.0129;
0083
0084   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 2;
0085   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 1;
0086   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 1;
0087   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 1;
0088
0089   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 3;
0090   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 1;
0091   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 1;
0092   cuts_[EleIDCutNames::MISSINGHITS][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 1;
0093
0094   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::VETO][EleIDEtaBins::BARREL] = 1;
0095   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::LOOSE][EleIDEtaBins::BARREL] = 1;
0096   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::BARREL] = 1;
0097   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::TIGHT][EleIDEtaBins::BARREL] = 1;
0098
0099   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::VETO][EleIDEtaBins::ENDCAP] = 1;
0100   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::LOOSE][EleIDEtaBins::ENDCAP] = 1;
0101   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::MEDIUM][EleIDEtaBins::ENDCAP] = 1;
0102   cuts_[EleIDCutNames::CONVERSION][EleIDWorkingPoints::TIGHT][EleIDEtaBins::ENDCAP] = 1;
0103 }
0104
0105 void ElectronIdentifier::setRho(double rho) {
0106   if (rho >= 0) {
0107     rho_ = rho;
0108   } else {
0109     throw cms::Exception("ValueError") << "Encountered invalid value for energy density rho.\n"
0110                                        << "Value: " << rho << "\n"
0111                                        << "Rho should be a real, positive number.\n";
0112   }
0113 }
0114 void ElectronIdentifier::setID(std::string ID) {
0115   if (ID == "TIGHT")
0116     ID_ = EleIDWorkingPoints::TIGHT;
0117   else if (ID == "MEDIUM")
0118     ID_ = EleIDWorkingPoints::MEDIUM;
0119   else if (ID == "LOOSE")
0120     ID_ = EleIDWorkingPoints::LOOSE;
0121   else if (ID == "VETO")
0122     ID_ = EleIDWorkingPoints::VETO;
0123   else
0124     throw;
0125 }
0126 float ElectronIdentifier::dEtaInSeed(const reco::GsfElectronPtr& ele) {
0127   return ele->superCluster().isNonnull() && ele->superCluster()->seed().isNonnull()
0128              ? ele->deltaEtaSuperClusterTrackAtVtx() - ele->superCluster()->eta() + ele->superCluster()->seed()->eta()
0129              : std::numeric_limits<float>::max();
0130 }
0131 float ElectronIdentifier::isolation(const reco::GsfElectronPtr& ele) {
0132   if (rho_ < 0) {
0133     throw;
0134   }
0135   const reco::GsfElectron::PflowIsolationVariables& pfIso = ele->pfIsolationVariables();
0136   const float chad = pfIso.sumChargedHadronPt;
0137   const float nhad = pfIso.sumNeutralHadronEt;
0138   const float pho = pfIso.sumPhotonEt;
0139   const float eA = _effectiveAreas.getEffectiveArea(fabs(ele->superCluster()->eta()));
0140   const float iso = chad + std::max(0.0, nhad + pho - rho_ * eA);
0141
0142   // Apply the cut and return the result
0143   // Scale by pT if the relative isolation is requested but avoid division by 0
0144   return iso;
0145 }
0146
0147 bool ElectronIdentifier::passID(const reco::GsfElectronPtr& ele,
0148                                 edm::Handle<reco::BeamSpot> beamspot,
0149                                 edm::Handle<reco::ConversionCollection> conversions) {
0150   if (ID_ == -1)
0151     throw;
0152   unsigned int region = fabs(ele->superCluster()->eta()) < 1.479 ? EleIDEtaBins::BARREL : EleIDEtaBins::BARREL;
0153
0154   if (ele->full5x5_sigmaIetaIeta() > cuts_[EleIDCutNames::SIGMAIETA][ID_][region])
0155     return false;
0156   if (dEtaInSeed(ele) > cuts_[EleIDCutNames::DETAINSEED][ID_][region])
0157     return false;
0158   if (std::abs(ele->deltaPhiSuperClusterTrackAtVtx()) > cuts_[EleIDCutNames::DPHIIN][ID_][region])
0159     return false;
0160   if (ele->hadronicOverEm() > cuts_[EleIDCutNames::HOVERE][ID_][region])
0161     return false;
0162   if (isolation(ele) / ele->pt() > cuts_[EleIDCutNames::ISO][ID_][region])
0163     return false;
0164   if (std::abs(1.0 - ele->eSuperClusterOverP()) / ele->ecalEnergy() > cuts_[EleIDCutNames::ONEOVERE][ID_][region])
0165     return false;
0166   if ((ele->gsfTrack()->hitPattern().numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS)) >
0167       cuts_[EleIDCutNames::MISSINGHITS][ID_][region])
0168     return false;
0169   if (ConversionTools::hasMatchedConversion(*ele, *conversions, beamspot->position()))
0170     return false;
0171
0172   return true;
0173 }