ElectronTools/src/EGammaCutBasedEleId.cc

0001 #include "EgammaAnalysis/ElectronTools/interface/EGammaCutBasedEleId.h"
0002 #include "EgammaAnalysis/ElectronTools/interface/ElectronEffectiveArea.h"
0003 #include "CommonTools/Egamma/interface/ConversionTools.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005
0006 #include <algorithm>
0007
0008 #ifndef STANDALONEID
0009
0010 bool EgammaCutBasedEleId::PassWP(WorkingPoint workingPoint,
0011                                  const reco::GsfElectron &ele,
0012                                  const edm::Handle<reco::ConversionCollection> &conversions,
0013                                  const reco::BeamSpot &beamspot,
0014                                  const edm::Handle<reco::VertexCollection> &vtxs,
0015                                  const double &iso_ch,
0016                                  const double &iso_em,
0017                                  const double &iso_nh,
0018                                  const double &rho,
0019                                  ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0020   // get the mask
0021   unsigned int mask = TestWP(workingPoint, ele, conversions, beamspot, vtxs, iso_ch, iso_em, iso_nh, rho, EAtarget);
0022
0023   // check if the desired WP passed
0024   if ((mask & PassAll) == PassAll)
0025     return true;
0026   return false;
0027 }
0028
0029 bool EgammaCutBasedEleId::PassWP(WorkingPoint workingPoint,
0030                                  const reco::GsfElectronRef &ele,
0031                                  const edm::Handle<reco::ConversionCollection> &conversions,
0032                                  const reco::BeamSpot &beamspot,
0033                                  const edm::Handle<reco::VertexCollection> &vtxs,
0034                                  const double &iso_ch,
0035                                  const double &iso_em,
0036                                  const double &iso_nh,
0037                                  const double &rho,
0038                                  ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0039   return PassWP(workingPoint, *ele, conversions, beamspot, vtxs, iso_ch, iso_em, iso_nh, rho, EAtarget);
0040 }
0041
0042 bool EgammaCutBasedEleId::PassTriggerCuts(TriggerWorkingPoint triggerWorkingPoint, const reco::GsfElectron &ele) {
0043   // get the variables
0044   bool isEB = ele.isEB() ? true : false;
0045   float pt = ele.pt();
0046   float dEtaIn = ele.deltaEtaSuperClusterTrackAtVtx();
0047   float dPhiIn = ele.deltaPhiSuperClusterTrackAtVtx();
0048   float sigmaIEtaIEta = ele.sigmaIetaIeta();
0049   float hoe = ele.hadronicOverEm();
0050   float trackIso = ele.dr03TkSumPt();
0051   float ecalIso = ele.dr03EcalRecHitSumEt();
0052   float hcalIso = ele.dr03HcalTowerSumEt();
0053
0054   // test the trigger cuts
0055   return EgammaCutBasedEleId::PassTriggerCuts(
0056       triggerWorkingPoint, isEB, pt, dEtaIn, dPhiIn, sigmaIEtaIEta, hoe, trackIso, ecalIso, hcalIso);
0057 }
0058
0059 bool EgammaCutBasedEleId::PassTriggerCuts(TriggerWorkingPoint triggerWorkingPoint, const reco::GsfElectronRef &ele) {
0060   return EgammaCutBasedEleId::PassTriggerCuts(triggerWorkingPoint, *ele);
0061 }
0062
0063 bool EgammaCutBasedEleId::PassEoverPCuts(const reco::GsfElectron &ele) {
0064   // get the variables
0065   float eta = ele.superCluster()->eta();
0066   float eopin = ele.eSuperClusterOverP();
0067   float fbrem = ele.fbrem();
0068
0069   // test the eop/fbrem cuts
0070   return EgammaCutBasedEleId::PassEoverPCuts(eta, eopin, fbrem);
0071 }
0072
0073 bool EgammaCutBasedEleId::PassEoverPCuts(const reco::GsfElectronRef &ele) { return PassEoverPCuts(*ele); }
0074
0075 unsigned int EgammaCutBasedEleId::TestWP(WorkingPoint workingPoint,
0076                                          const reco::GsfElectron &ele,
0077                                          const edm::Handle<reco::ConversionCollection> &conversions,
0078                                          const reco::BeamSpot &beamspot,
0079                                          const edm::Handle<reco::VertexCollection> &vtxs,
0080                                          const double &iso_ch,
0081                                          const double &iso_em,
0082                                          const double &iso_nh,
0083                                          const double &rho,
0084                                          ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0085   // get the ID variables from the electron object
0086
0087   // kinematic variables
0088   bool isEB = ele.isEB() ? true : false;
0089   float pt = ele.pt();
0090   float eta = ele.superCluster()->eta();
0091
0092   // id variables
0093   float dEtaIn = ele.deltaEtaSuperClusterTrackAtVtx();
0094   float dPhiIn = ele.deltaPhiSuperClusterTrackAtVtx();
0095   float sigmaIEtaIEta = ele.sigmaIetaIeta();
0096   float hoe = ele.hadronicOverEm();
0097   float ooemoop = (1.0 / ele.ecalEnergy() - ele.eSuperClusterOverP() / ele.ecalEnergy());
0098
0099   // impact parameter variables
0100   float d0vtx = 0.0;
0101   float dzvtx = 0.0;
0102   if (!vtxs->empty()) {
0103     reco::VertexRef vtx(vtxs, 0);
0104     d0vtx = ele.gsfTrack()->dxy(vtx->position());
0105     dzvtx = ele.gsfTrack()->dz(vtx->position());
0106   } else {
0107     d0vtx = ele.gsfTrack()->dxy();
0108     dzvtx = ele.gsfTrack()->dz();
0109   }
0110
0111   // conversion rejection variables
0112   bool vtxFitConversion = ConversionTools::hasMatchedConversion(ele, *conversions, beamspot.position());
0113   float mHits = ele.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0114
0115   // get the mask value
0116   unsigned int mask = EgammaCutBasedEleId::TestWP(workingPoint,
0117                                                   isEB,
0118                                                   pt,
0119                                                   eta,
0120                                                   dEtaIn,
0121                                                   dPhiIn,
0122                                                   sigmaIEtaIEta,
0123                                                   hoe,
0124                                                   ooemoop,
0125                                                   d0vtx,
0126                                                   dzvtx,
0127                                                   iso_ch,
0128                                                   iso_em,
0129                                                   iso_nh,
0130                                                   vtxFitConversion,
0131                                                   mHits,
0132                                                   rho,
0133                                                   EAtarget);
0134
0135   // return the mask value
0136   return mask;
0137 }
0138
0139 unsigned int EgammaCutBasedEleId::TestWP(WorkingPoint workingPoint,
0140                                          const reco::GsfElectronRef &ele,
0141                                          const edm::Handle<reco::ConversionCollection> &conversions,
0142                                          const reco::BeamSpot &beamspot,
0143                                          const edm::Handle<reco::VertexCollection> &vtxs,
0144                                          const double &iso_ch,
0145                                          const double &iso_em,
0146                                          const double &iso_nh,
0147                                          const double &rho,
0148                                          ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0149   return TestWP(workingPoint, *ele, conversions, beamspot, vtxs, iso_ch, iso_em, iso_nh, rho, EAtarget);
0150 }
0151
0152 #endif
0153
0154 bool EgammaCutBasedEleId::PassWP(WorkingPoint workingPoint,
0155                                  const bool isEB,
0156                                  const float pt,
0157                                  const float eta,
0158                                  const float dEtaIn,
0159                                  const float dPhiIn,
0160                                  const float sigmaIEtaIEta,
0161                                  const float hoe,
0162                                  const float ooemoop,
0163                                  const float d0vtx,
0164                                  const float dzvtx,
0165                                  const float iso_ch,
0166                                  const float iso_em,
0167                                  const float iso_nh,
0168                                  const bool vtxFitConversion,
0169                                  const unsigned int mHits,
0170                                  const double rho,
0171                                  ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0172   unsigned int mask = EgammaCutBasedEleId::TestWP(workingPoint,
0173                                                   isEB,
0174                                                   pt,
0175                                                   eta,
0176                                                   dEtaIn,
0177                                                   dPhiIn,
0178                                                   sigmaIEtaIEta,
0179                                                   hoe,
0180                                                   ooemoop,
0181                                                   d0vtx,
0182                                                   dzvtx,
0183                                                   iso_ch,
0184                                                   iso_em,
0185                                                   iso_nh,
0186                                                   vtxFitConversion,
0187                                                   mHits,
0188                                                   rho,
0189                                                   EAtarget);
0190
0191   if ((mask & PassAll) == PassAll)
0192     return true;
0193   return false;
0194 }
0195
0196 bool EgammaCutBasedEleId::PassTriggerCuts(const TriggerWorkingPoint triggerWorkingPoint,
0197                                           const bool isEB,
0198                                           const float pt,
0199                                           const float dEtaIn,
0200                                           const float dPhiIn,
0201                                           const float sigmaIEtaIEta,
0202                                           const float hoe,
0203                                           const float trackIso,
0204                                           const float ecalIso,
0205                                           const float hcalIso) {
0206   // choose cut if barrel or endcap
0207   unsigned int idx = isEB ? 0 : 1;
0208
0209   if (triggerWorkingPoint == EgammaCutBasedEleId::TRIGGERTIGHT) {
0210     float cut_dEtaIn[2] = {0.007, 0.009};
0211     float cut_dPhiIn[2] = {0.15, 0.10};
0212     float cut_sigmaIEtaIEta[2] = {0.01, 0.03};
0213     float cut_hoe[2] = {0.12, 0.10};
0214     float cut_trackIso[2] = {0.20, 0.20};
0215     float cut_ecalIso[2] = {0.20, 0.20};
0216     float cut_hcalIso[2] = {0.20, 0.20};
0217     if (fabs(dEtaIn) > cut_dEtaIn[idx])
0218       return false;
0219     if (fabs(dPhiIn) > cut_dPhiIn[idx])
0220       return false;
0221     if (sigmaIEtaIEta > cut_sigmaIEtaIEta[idx])
0222       return false;
0223     if (hoe > cut_hoe[idx])
0224       return false;
0225     if (trackIso / pt > cut_trackIso[idx])
0226       return false;
0227     if (ecalIso / pt > cut_ecalIso[idx])
0228       return false;
0229     if (hcalIso / pt > cut_hcalIso[idx])
0230       return false;
0231   } else if (triggerWorkingPoint == EgammaCutBasedEleId::TRIGGERWP70) {
0232     float cut_dEtaIn[2] = {0.004, 0.005};
0233     float cut_dPhiIn[2] = {0.03, 0.02};
0234     float cut_sigmaIEtaIEta[2] = {0.01, 0.03};
0235     float cut_hoe[2] = {0.025, 0.025};
0236     float cut_trackIso[2] = {0.10, 0.10};
0237     float cut_ecalIso[2] = {0.10, 0.05};
0238     float cut_hcalIso[2] = {0.05, 0.05};
0239     if (fabs(dEtaIn) > cut_dEtaIn[idx])
0240       return false;
0241     if (fabs(dPhiIn) > cut_dPhiIn[idx])
0242       return false;
0243     if (sigmaIEtaIEta > cut_sigmaIEtaIEta[idx])
0244       return false;
0245     if (hoe > cut_hoe[idx])
0246       return false;
0247     if (trackIso / pt > cut_trackIso[idx])
0248       return false;
0249     if (ecalIso / pt > cut_ecalIso[idx])
0250       return false;
0251     if (hcalIso / pt > cut_hcalIso[idx])
0252       return false;
0253   } else {
0254     std::cout << "[EgammaCutBasedEleId::PassTriggerCuts] Undefined working point" << std::endl;
0255   }
0256
0257   return true;
0258 }
0259
0260 bool EgammaCutBasedEleId::PassEoverPCuts(const float eta, const float eopin, const float fbrem) {
0261   if (fbrem > 0.15)
0262     return true;
0263   else if (fabs(eta) < 1.0 && eopin > 0.95)
0264     return true;
0265   return false;
0266 }
0267
0268 unsigned int EgammaCutBasedEleId::TestWP(WorkingPoint workingPoint,
0269                                          const bool isEB,
0270                                          const float pt,
0271                                          const float eta,
0272                                          const float dEtaIn,
0273                                          const float dPhiIn,
0274                                          const float sigmaIEtaIEta,
0275                                          const float hoe,
0276                                          const float ooemoop,
0277                                          const float d0vtx,
0278                                          const float dzvtx,
0279                                          const float iso_ch,
0280                                          const float iso_em,
0281                                          const float iso_nh,
0282                                          const bool vtxFitConversion,
0283                                          const unsigned int mHits,
0284                                          const double rho,
0285                                          ElectronEffectiveArea::ElectronEffectiveAreaTarget EAtarget) {
0286   unsigned int mask = 0;
0287   float cut_dEtaIn[2] = {999.9, 999.9};
0288   float cut_dPhiIn[2] = {999.9, 999.9};
0289   float cut_sigmaIEtaIEta[2] = {999.9, 999.9};
0290   float cut_hoe[2] = {999.9, 999.9};
0291   float cut_ooemoop[2] = {999.9, 999.9};
0292   float cut_d0vtx[2] = {999.9, 999.9};
0293   float cut_dzvtx[2] = {999.9, 999.9};
0294   float cut_iso[2] = {999.9, 999.9};
0295   bool cut_vtxFit[2] = {false, false};
0296   unsigned int cut_mHits[2] = {999, 999};
0297
0298   if (workingPoint == EgammaCutBasedEleId::VETO) {
0299     cut_dEtaIn[0] = 0.007;
0300     cut_dEtaIn[1] = 0.010;
0301     cut_dPhiIn[0] = 0.800;
0302     cut_dPhiIn[1] = 0.700;
0303     cut_sigmaIEtaIEta[0] = 0.010;
0304     cut_sigmaIEtaIEta[1] = 0.030;
0305     cut_hoe[0] = 0.150;
0306     cut_hoe[1] = 999.9;
0307     cut_ooemoop[0] = 999.9;
0308     cut_ooemoop[1] = 999.9;
0309     cut_d0vtx[0] = 0.040;
0310     cut_d0vtx[1] = 0.040;
0311     cut_dzvtx[0] = 0.200;
0312     cut_dzvtx[1] = 0.200;
0313     cut_vtxFit[0] = false;
0314     cut_vtxFit[1] = false;
0315     cut_mHits[0] = 999;
0316     cut_mHits[1] = 999;
0317     cut_iso[0] = 0.150;
0318     cut_iso[1] = 0.150;
0319   } else if (workingPoint == EgammaCutBasedEleId::LOOSE) {
0320     cut_dEtaIn[0] = 0.007;
0321     cut_dEtaIn[1] = 0.009;
0322     cut_dPhiIn[0] = 0.150;
0323     cut_dPhiIn[1] = 0.100;
0324     cut_sigmaIEtaIEta[0] = 0.010;
0325     cut_sigmaIEtaIEta[1] = 0.030;
0326     cut_hoe[0] = 0.120;
0327     cut_hoe[1] = 0.100;
0328     cut_ooemoop[0] = 0.050;
0329     cut_ooemoop[1] = 0.050;
0330     cut_d0vtx[0] = 0.020;
0331     cut_d0vtx[1] = 0.020;
0332     cut_dzvtx[0] = 0.200;
0333     cut_dzvtx[1] = 0.200;
0334     cut_vtxFit[0] = true;
0335     cut_vtxFit[1] = true;
0336     cut_mHits[0] = 1;
0337     cut_mHits[1] = 1;
0338     if (pt >= 20.0) {
0339       cut_iso[0] = 0.150;
0340       cut_iso[1] = 0.150;
0341     } else {
0342       cut_iso[0] = 0.150;
0343       cut_iso[1] = 0.100;
0344     }
0345   } else if (workingPoint == EgammaCutBasedEleId::MEDIUM) {
0346     cut_dEtaIn[0] = 0.004;
0347     cut_dEtaIn[1] = 0.007;
0348     cut_dPhiIn[0] = 0.060;
0349     cut_dPhiIn[1] = 0.030;
0350     cut_sigmaIEtaIEta[0] = 0.010;
0351     cut_sigmaIEtaIEta[1] = 0.030;
0352     cut_hoe[0] = 0.120;
0353     cut_hoe[1] = 0.100;
0354     cut_ooemoop[0] = 0.050;
0355     cut_ooemoop[1] = 0.050;
0356     cut_d0vtx[0] = 0.020;
0357     cut_d0vtx[1] = 0.020;
0358     cut_dzvtx[0] = 0.100;
0359     cut_dzvtx[1] = 0.100;
0360     cut_vtxFit[0] = true;
0361     cut_vtxFit[1] = true;
0362     cut_mHits[0] = 1;
0363     cut_mHits[1] = 1;
0364     if (pt >= 20.0) {
0365       cut_iso[0] = 0.150;
0366       cut_iso[1] = 0.150;
0367     } else {
0368       cut_iso[0] = 0.150;
0369       cut_iso[1] = 0.100;
0370     }
0371   } else if (workingPoint == EgammaCutBasedEleId::TIGHT) {
0372     cut_dEtaIn[0] = 0.004;
0373     cut_dEtaIn[1] = 0.005;
0374     cut_dPhiIn[0] = 0.030;
0375     cut_dPhiIn[1] = 0.020;
0376     cut_sigmaIEtaIEta[0] = 0.010;
0377     cut_sigmaIEtaIEta[1] = 0.030;
0378     cut_hoe[0] = 0.120;
0379     cut_hoe[1] = 0.100;
0380     cut_ooemoop[0] = 0.050;
0381     cut_ooemoop[1] = 0.050;
0382     cut_d0vtx[0] = 0.020;
0383     cut_d0vtx[1] = 0.020;
0384     cut_dzvtx[0] = 0.100;
0385     cut_dzvtx[1] = 0.100;
0386     cut_vtxFit[0] = true;
0387     cut_vtxFit[1] = true;
0388     cut_mHits[0] = 0;
0389     cut_mHits[1] = 0;
0390     if (pt >= 20.0) {
0391       cut_iso[0] = 0.100;
0392       cut_iso[1] = 0.100;
0393     } else {
0394       cut_iso[0] = 0.100;
0395       cut_iso[1] = 0.070;
0396     }
0397   } else {
0398     std::cout << "[EgammaCutBasedEleId::TestWP] Undefined working point" << std::endl;
0399   }
0400
0401   // choose cut if barrel or endcap
0402   unsigned int idx = isEB ? 0 : 1;
0403
0404   // effective area for isolation
0405   float AEff = ElectronEffectiveArea::GetElectronEffectiveArea(
0406       ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, eta, EAtarget);
0407   //float AEff = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, eta, ElectronEffectiveArea::kEleEAData2011);
0408
0409   // apply to neutrals
0410   double rhoPrime = std::max(rho, 0.0);
0411   double iso_n = std::max(iso_nh + iso_em - rhoPrime * AEff, 0.0);
0412
0413   // compute final isolation
0414   double iso = (iso_n + iso_ch) / pt;
0415
0416   // test cuts
0417   if (fabs(dEtaIn) < cut_dEtaIn[idx])
0418     mask |= DETAIN;
0419   if (fabs(dPhiIn) < cut_dPhiIn[idx])
0420     mask |= DPHIIN;
0421   if (sigmaIEtaIEta < cut_sigmaIEtaIEta[idx])
0422     mask |= SIGMAIETAIETA;
0423   if (hoe < cut_hoe[idx])
0424     mask |= HOE;
0425   if (fabs(ooemoop) < cut_ooemoop[idx])
0426     mask |= OOEMOOP;
0427   if (fabs(d0vtx) < cut_d0vtx[idx])
0428     mask |= D0VTX;
0429   if (fabs(dzvtx) < cut_dzvtx[idx])
0430     mask |= DZVTX;
0431   if (!cut_vtxFit[idx] || !vtxFitConversion)
0432     mask |= VTXFIT;
0433   if (mHits <= cut_mHits[idx])
0434     mask |= MHITS;
0435   if (iso < cut_iso[idx])
0436     mask |= ISO;
0437
0438   // return the mask
0439   return mask;
0440 }
0441
0442 void EgammaCutBasedEleId::PrintDebug(unsigned int mask) {
0443   printf("detain(%i), ", bool(mask & DETAIN));
0444   printf("dphiin(%i), ", bool(mask & DPHIIN));
0445   printf("sieie(%i), ", bool(mask & SIGMAIETAIETA));
0446   printf("hoe(%i), ", bool(mask & HOE));
0447   printf("ooemoop(%i), ", bool(mask & OOEMOOP));
0448   printf("d0vtx(%i), ", bool(mask & D0VTX));
0449   printf("dzvtx(%i), ", bool(mask & DZVTX));
0450   printf("iso(%i), ", bool(mask & ISO));
0451   printf("vtxfit(%i), ", bool(mask & VTXFIT));
0452   printf("mhits(%i)\n", bool(mask & MHITS));
0453 }