LowPtElectrons: updated BDT model and code base for ID

RecoEgamma/EgammaElectronProducers/BuildFile.xml

@@ -0,0 +1,11 @@
+<use name="CommonTools/BaseParticlePropagator"/>
+<use name="DataFormats/BeamSpot"/>
+<use name="DataFormats/Common"/>
+<use name="DataFormats/EgammaCandidates"/>
+<use name="DataFormats/GsfTrackReco"/>
+<use name="DataFormats/ParticleFlowReco"/>
+<use name="DataFormats/TrackReco"/>
+<use name="RecoEcal/EgammaCoreTools"/>
+<export>
+  <lib name="RecoEgammaEgammaElectronProducers"/>
+</export>

RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h

@@ -0,0 +1,34 @@
+#ifndef RecoEgamma_EgammaElectronProducers_LowPtGsfElectronFeatures_h
+#define RecoEgamma_EgammaElectronProducers_LowPtGsfElectronFeatures_h
+
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
+#include "DataFormats/Common/interface/Ptr.h"
+#include "DataFormats/Common/interface/Ref.h"
+#include "DataFormats/Common/interface/RefToPtr.h"
+#include "DataFormats/Common/interface/View.h"
+#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
+#include "DataFormats/ParticleFlowReco/interface/PreId.h"
+#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
+#include <vector>
+
+namespace lowptgsfeleseed {
+
+  std::vector<float> features(const reco::PreId& ecal,
+                              const reco::PreId& hcal,
+                              double rho,
+                              const reco::BeamSpot& spot,
+                              noZS::EcalClusterLazyTools& ecalTools);
+
+}
+
+namespace lowptgsfeleid {
+
+  std::vector<float> features(edm::Ptr<reco::GsfElectron> const& ele, float rho, float unbiased);
+
+  std::vector<float> features(edm::Ref<std::vector<reco::GsfElectron> > const& ele, float rho, float unbiased);
+
+  std::vector<float> features(edm::Ref<edm::View<reco::GsfElectron> > const& ele, float rho, float unbiased);
+
+}  // namespace lowptgsfeleid
+
+#endif  // RecoEgamma_EgammaElectronProducers_LowPtGsfElectronFeatures_h

RecoEgamma/EgammaElectronProducers/plugins/BuildFile.xml

@@ -20,6 +20,7 @@
 <use name="MagneticField/Records"/>
 <use name="RecoEcal/EgammaCoreTools"/>
 <use name="RecoEgamma/EgammaElectronAlgos"/>
+<use name="RecoEgamma/EgammaElectronProducers"/>
 <use name="RecoEgamma/EgammaIsolationAlgos"/>
 <use name="RecoParticleFlow/PFClusterTools"/>
 <use name="RecoTracker/TkTrackingRegions"/>

RecoEgamma/EgammaElectronProducers/plugins/LowPtGsfElectronIDProducer.cc

@@ -1,5 +1,8 @@
+#include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
 #include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/Common/interface/Ptr.h"
 #include "DataFormats/Common/interface/ValueMap.h"
+#include "DataFormats/Common/interface/View.h"
 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
 #include "FWCore/Utilities/interface/InputTag.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
@@ -11,124 +14,14 @@
 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
 #include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
 #include "FWCore/Framework/interface/Event.h"
-
+#include "RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h"
 #include <string>
 #include <vector>
 
-namespace {
-
-  std::vector<float> getFeatures(reco::GsfElectronRef const& ele, float rho) {
-    // KF track
-    float trk_p = -1.;
-    float trk_nhits = -1.;
-    float trk_chi2red = -1.;
-    // GSF track
-    float gsf_nhits = -1.;
-    float gsf_chi2red = -1.;
-    // SC
-    float sc_E = -1.;
-    float sc_eta = -1.;
-    float sc_etaWidth = -1.;
-    float sc_phiWidth = -1.;
-    // Track-cluster matching
-    float match_seed_dEta = -1.;
-    float match_eclu_EoverP = -1.;
-    float match_SC_EoverP = -1.;
-    float match_SC_dEta = -1.;
-    float match_SC_dPhi = -1.;
-    // Shower shape vars
-    float shape_full5x5_sigmaIetaIeta = -1.;
-    float shape_full5x5_sigmaIphiIphi = -1.;
-    float shape_full5x5_HoverE = -1.;
-    float shape_full5x5_r9 = -1.;
-    float shape_full5x5_circularity = -1.;
-    // Misc
-    float brem_frac = -1.;
-    float ele_pt = -1.;
-
-    // KF tracks
-    if (ele->core().isNonnull()) {
-      reco::TrackRef trk = ele->core()->ctfTrack();  //@@ is this what we want?!
-      if (trk.isNonnull()) {
-        trk_p = float(trk->p());
-        trk_nhits = float(trk->found());
-        trk_chi2red = float(trk->normalizedChi2());
-      }
-    }
-
-    // GSF tracks
-    if (ele->core().isNonnull()) {
-      reco::GsfTrackRef gsf = ele->core()->gsfTrack();
-      if (gsf.isNonnull()) {
-        gsf_nhits = gsf->found();
-        gsf_chi2red = gsf->normalizedChi2();
-      }
-    }
-
-    // Super clusters
-    if (ele->core().isNonnull()) {
-      reco::SuperClusterRef sc = ele->core()->superCluster();
-      if (sc.isNonnull()) {
-        sc_E = sc->energy();
-        sc_eta = sc->eta();
-        sc_etaWidth = sc->etaWidth();
-        sc_phiWidth = sc->phiWidth();
-      }
-    }
-
-    // Track-cluster matching
-    if (ele.isNonnull()) {
-      match_seed_dEta = ele->deltaEtaSeedClusterTrackAtCalo();
-      match_eclu_EoverP = (1. / ele->ecalEnergy()) - (1. / ele->p());
-      match_SC_EoverP = ele->eSuperClusterOverP();
-      match_SC_dEta = ele->deltaEtaSuperClusterTrackAtVtx();
-      match_SC_dPhi = ele->deltaPhiSuperClusterTrackAtVtx();
-    }
-
-    // Shower shape vars
-    if (ele.isNonnull()) {
-      shape_full5x5_sigmaIetaIeta = ele->full5x5_sigmaIetaIeta();
-      shape_full5x5_sigmaIphiIphi = ele->full5x5_sigmaIphiIphi();
-      shape_full5x5_HoverE = ele->full5x5_hcalOverEcal();
-      shape_full5x5_r9 = ele->full5x5_r9();
-      shape_full5x5_circularity = 1. - ele->full5x5_e1x5() / ele->full5x5_e5x5();
-    }
-
-    // Misc
-    if (ele.isNonnull()) {
-      brem_frac = ele->fbrem();
-      ele_pt = ele->pt();
-    }
-
-    return {
-        rho,
-        ele_pt,
-        sc_eta,
-        shape_full5x5_sigmaIetaIeta,
-        shape_full5x5_sigmaIphiIphi,
-        shape_full5x5_circularity,
-        shape_full5x5_r9,
-        sc_etaWidth,
-        sc_phiWidth,
-        shape_full5x5_HoverE,
-        trk_nhits,
-        trk_chi2red,
-        gsf_chi2red,
-        brem_frac,
-        gsf_nhits,
-        match_SC_EoverP,
-        match_eclu_EoverP,
-        match_SC_dEta,
-        match_SC_dPhi,
-        match_seed_dEta,
-        sc_E,
-        trk_p,
-    };
-  }
-
-}  // namespace
-
+////////////////////////////////////////////////////////////////////////////////
+//
 class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> {
 public:
   explicit LowPtGsfElectronIDProducer(const edm::ParameterSet&);
@@ -138,10 +31,11 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> {
   static void fillDescriptions(edm::ConfigurationDescriptions&);
 
 private:
-  double eval(const std::string& name, const reco::GsfElectronRef&, double rho) const;
+  double eval(const std::string& name, const edm::Ptr<reco::GsfElectron>&, double rho, float unbiased) const;
 
-  const edm::EDGetTokenT<reco::GsfElectronCollection> gsfElectrons_;
+  const edm::EDGetTokenT<edm::View<reco::GsfElectron> > gsfElectrons_;
   const edm::EDGetTokenT<double> rho_;
+  const edm::EDGetTokenT<edm::ValueMap<float> > unbiased_;
   const std::vector<std::string> names_;
   const bool passThrough_;
   const double minPtThreshold_;
@@ -153,8 +47,9 @@ class LowPtGsfElectronIDProducer final : public edm::global::EDProducer<> {
 ////////////////////////////////////////////////////////////////////////////////
 //
 LowPtGsfElectronIDProducer::LowPtGsfElectronIDProducer(const edm::ParameterSet& conf)
-    : gsfElectrons_(consumes<reco::GsfElectronCollection>(conf.getParameter<edm::InputTag>("electrons"))),
+    : gsfElectrons_(consumes<edm::View<reco::GsfElectron> >(conf.getParameter<edm::InputTag>("electrons"))),
       rho_(consumes<double>(conf.getParameter<edm::InputTag>("rho"))),
+      unbiased_(consumes<edm::ValueMap<float> >(conf.getParameter<edm::InputTag>("unbiased"))),
       names_(conf.getParameter<std::vector<std::string> >("ModelNames")),
       passThrough_(conf.getParameter<bool>("PassThrough")),
       minPtThreshold_(conf.getParameter<double>("MinPtThreshold")),
@@ -183,14 +78,30 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const
   edm::Handle<double> rho;
   event.getByToken(rho_, rho);
   if (!rho.isValid()) {
-    edm::LogError("Problem with rho handle");
+    std::ostringstream os;
+    os << "Problem accessing rho collection for low-pT electrons" << std::endl;
+    edm::LogError("InvalidHandle") << os.str();
+    throw cms::Exception("InvalidHandle", os.str());
   }
 
   // Retrieve GsfElectrons from Event
-  edm::Handle<reco::GsfElectronCollection> gsfElectrons;
+  edm::Handle<edm::View<reco::GsfElectron> > gsfElectrons;
   event.getByToken(gsfElectrons_, gsfElectrons);
   if (!gsfElectrons.isValid()) {
-    edm::LogError("Problem with gsfElectrons handle");
+    std::ostringstream os;
+    os << "Problem accessing low-pT gsfElectrons collection" << std::endl;
+    edm::LogError("InvalidHandle") << os.str();
+    throw cms::Exception("InvalidHandle", os.str());
+  }
+
+  // ElectronSeed unbiased BDT
+  edm::Handle<edm::ValueMap<float> > unbiasedH;
+  event.getByToken(unbiased_, unbiasedH);
+  if (!unbiasedH.isValid()) {
+    std::ostringstream os;
+    os << "Problem accessing low-pT 'unbiased' ElectronSeed collection" << std::endl;
+    edm::LogError("InvalidHandle") << os.str();
+    throw cms::Exception("InvalidHandle", os.str());
   }
 
   // Iterate through Electrons, evaluate BDT, and store result
@@ -199,10 +110,20 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const
     output.emplace_back(gsfElectrons->size(), -999.);
   }
   for (unsigned int iele = 0; iele < gsfElectrons->size(); iele++) {
-    reco::GsfElectronRef ele(gsfElectrons, iele);
+    edm::Ptr<reco::GsfElectron> ele(gsfElectrons, iele);
+
+    if (ele->core().isNull()) {
+      continue;
+    }
+    const auto& gsf = ele->core()->gsfTrack();  // reco::GsfTrackRef
+    if (gsf.isNull()) {
+      continue;
+    }
+    float unbiased = (*unbiasedH)[gsf];
+
     //if ( !passThrough_ && ( ele->pt() < minPtThreshold_ ) ) { continue; }
     for (unsigned int iname = 0; iname < names_.size(); ++iname) {
-      output[iname][iele] = eval(names_[iname], ele, *rho);
+      output[iname][iele] = eval(names_[iname], ele, *rho, unbiased);
     }
   }
 
@@ -212,16 +133,20 @@ void LowPtGsfElectronIDProducer::produce(edm::StreamID, edm::Event& event, const
     edm::ValueMap<float>::Filler filler(*ptr);
     filler.insert(gsfElectrons, output[iname].begin(), output[iname].end());
     filler.fill();
-    reco::GsfElectronRef ele(gsfElectrons, 0);
     event.put(std::move(ptr), names_[iname]);
   }
 }
 
-double LowPtGsfElectronIDProducer::eval(const std::string& name, const reco::GsfElectronRef& ele, double rho) const {
+//////////////////////////////////////////////////////////////////////////////////////////
+//
+double LowPtGsfElectronIDProducer::eval(const std::string& name,
+                                        const edm::Ptr<reco::GsfElectron>& ele,
+                                        double rho,
+                                        float unbiased) const {
   auto iter = std::find(names_.begin(), names_.end(), name);
   if (iter != names_.end()) {
     int index = std::distance(names_.begin(), iter);
-    std::vector<float> inputs = getFeatures(ele, rho);
+    std::vector<float> inputs = lowptgsfeleid::features(ele, rho, unbiased);
     return models_.at(index)->GetResponse(inputs.data());
   } else {
     throw cms::Exception("Unknown model name") << "'Name given: '" << name << "'. Check against configuration file.\n";
@@ -234,6 +159,7 @@ double LowPtGsfElectronIDProducer::eval(const std::string& name, const reco::Gsf
 void LowPtGsfElectronIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
   desc.add<edm::InputTag>("electrons", edm::InputTag("lowPtGsfElectrons"));
+  desc.add<edm::InputTag>("unbiased", edm::InputTag("lowPtGsfElectronSeedValueMaps:unbiased"));
   desc.add<edm::InputTag>("rho", edm::InputTag("fixedGridRhoFastjetAllTmp"));
   desc.add<std::vector<std::string> >("ModelNames", {""});
   desc.add<std::vector<std::string> >(