HLT XGBoost Photon MVA and Diphoton combination filter (backport for 14_0_X)

HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

@@ -0,0 +1,157 @@
+#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "HLTrigger/HLTcore/interface/HLTFilter.h"
+
+#include <vector>
+#include <memory>
+
+namespace edm {
+  class ConfigurationDescriptions;
+}
+
+class HLTEgammaDoubleXGBoostCombFilter : public HLTFilter {
+public:
+  explicit HLTEgammaDoubleXGBoostCombFilter(edm::ParameterSet const&);
+  ~HLTEgammaDoubleXGBoostCombFilter() override {}
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  bool hltFilter(edm::Event& event,
+                 const edm::EventSetup& setup,
+                 trigger::TriggerFilterObjectWithRefs& filterproduct) const override;
+
+  const double highMassCut_;
+  const std::vector<double> leadCutHighMass1_;
+  const std::vector<double> subCutHighMass1_;
+  const std::vector<double> leadCutHighMass2_;
+  const std::vector<double> subCutHighMass2_;
+  const std::vector<double> leadCutHighMass3_;
+  const std::vector<double> subCutHighMass3_;
+
+  const double lowMassCut_;
+  const std::vector<double> leadCutLowMass1_;
+  const std::vector<double> subCutLowMass1_;
+  const std::vector<double> leadCutLowMass2_;
+  const std::vector<double> subCutLowMass2_;
+  const std::vector<double> leadCutLowMass3_;
+  const std::vector<double> subCutLowMass3_;
+
+  const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> candToken_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> mvaToken_;
+};
+
+HLTEgammaDoubleXGBoostCombFilter::HLTEgammaDoubleXGBoostCombFilter(edm::ParameterSet const& config)
+    : HLTFilter(config),
+      highMassCut_(config.getParameter<double>("highMassCut")),
+      leadCutHighMass1_(config.getParameter<std::vector<double>>("leadCutHighMass1")),
+      subCutHighMass1_(config.getParameter<std::vector<double>>("subCutHighMass1")),
+      leadCutHighMass2_(config.getParameter<std::vector<double>>("leadCutHighMass2")),
+      subCutHighMass2_(config.getParameter<std::vector<double>>("subCutHighMass2")),
+      leadCutHighMass3_(config.getParameter<std::vector<double>>("leadCutHighMass3")),
+      subCutHighMass3_(config.getParameter<std::vector<double>>("subCutHighMass3")),
+
+      lowMassCut_(config.getParameter<double>("lowMassCut")),
+      leadCutLowMass1_(config.getParameter<std::vector<double>>("leadCutLowMass1")),
+      subCutLowMass1_(config.getParameter<std::vector<double>>("subCutLowMass1")),
+      leadCutLowMass2_(config.getParameter<std::vector<double>>("leadCutLowMass2")),
+      subCutLowMass2_(config.getParameter<std::vector<double>>("subCutLowMass2")),
+      leadCutLowMass3_(config.getParameter<std::vector<double>>("leadCutLowMass3")),
+      subCutLowMass3_(config.getParameter<std::vector<double>>("subCutLowMass3")),
+
+      candToken_(consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("candTag"))),
+      mvaToken_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("mvaPhotonTag"))) {}
+
+void HLTEgammaDoubleXGBoostCombFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  makeHLTFilterDescription(desc);
+
+  desc.add<double>("highMassCut", 95.0);
+  desc.add<std::vector<double>>("leadCutHighMass1", {0.98, 0.95});
+  desc.add<std::vector<double>>("subCutHighMass1", {0.0, 0.04});
+  desc.add<std::vector<double>>("leadCutHighMass2", {0.85, 0.85});
+  desc.add<std::vector<double>>("subCutHighMass2", {0.04, 0.08});
+  desc.add<std::vector<double>>("leadCutHighMass3", {0.30, 0.50});
+  desc.add<std::vector<double>>("subCutHighMass3", {0.15, 0.20});
+
+  desc.add<double>("lowMassCut", 60.0);
+  desc.add<std::vector<double>>("leadCutLowMass1", {0.98, 0.90});
+  desc.add<std::vector<double>>("subCutLowMass1", {0.04, 0.05});
+  desc.add<std::vector<double>>("leadCutLowMass2", {0.90, 0.80});
+  desc.add<std::vector<double>>("subCutLowMass2", {0.10, 0.10});
+  desc.add<std::vector<double>>("leadCutLowMass3", {0.60, 0.60});
+  desc.add<std::vector<double>>("subCutLowMass3", {0.30, 0.30});
+
+  desc.add<edm::InputTag>("candTag", edm::InputTag("hltEgammaCandidatesUnseeded"));
+  desc.add<edm::InputTag>("mvaPhotonTag", edm::InputTag("PhotonXGBoostProducer"));
+
+  descriptions.addWithDefaultLabel(desc);
+}
+
+bool HLTEgammaDoubleXGBoostCombFilter::hltFilter(edm::Event& event,
+                                                 const edm::EventSetup& setup,
+                                                 trigger::TriggerFilterObjectWithRefs& filterproduct) const {
+  //producer collection (hltEgammaCandidates(Unseeded))
+  const auto& recCollection = event.getHandle(candToken_);
+
+  //get hold of photon MVA association map
+  const auto& mvaMap = event.getHandle(mvaToken_);
+
+  std::vector<math::XYZTLorentzVector> p4s(recCollection->size());
+  std::vector<bool> isTight(recCollection->size());
+
+  bool accept = false;
+
+  for (size_t i = 0; i < recCollection->size(); i++) {
+    edm::Ref<reco::RecoEcalCandidateCollection> refi(recCollection, i);
+    float EtaSCi = refi->eta();
+    int eta1 = (std::abs(EtaSCi) < 1.5) ? 0 : 1;
+    float mvaScorei = (*mvaMap).find(refi)->val;
+    math::XYZTLorentzVector p4i = refi->p4();
+    for (size_t j = i + 1; j < recCollection->size(); j++) {
+      edm::Ref<reco::RecoEcalCandidateCollection> refj(recCollection, j);
+      float EtaSCj = refj->eta();
+      int eta2 = (std::abs(EtaSCj) < 1.5) ? 0 : 1;
+      float mvaScorej = (*mvaMap).find(refj)->val;
+      math::XYZTLorentzVector p4j = refj->p4();
+      math::XYZTLorentzVector pairP4 = p4i + p4j;
+      double mass = pairP4.M();
+      if (mass >= highMassCut_) {
+        if (mvaScorei >= mvaScorej && ((mvaScorei > leadCutHighMass1_[eta1] && mvaScorej > subCutHighMass1_[eta2]) ||
+                                       (mvaScorei > leadCutHighMass2_[eta1] && mvaScorej > subCutHighMass2_[eta2]) ||
+                                       (mvaScorei > leadCutHighMass3_[eta1] && mvaScorej > subCutHighMass3_[eta2]))) {
+          accept = true;
+        }  //if scoreI > scoreJ
+        else if (mvaScorej > mvaScorei &&
+                 ((mvaScorej > leadCutHighMass1_[eta1] && mvaScorei > subCutHighMass1_[eta2]) ||
+                  (mvaScorej > leadCutHighMass2_[eta1] && mvaScorei > subCutHighMass2_[eta2]) ||
+                  (mvaScorej > leadCutHighMass3_[eta1] && mvaScorei > subCutHighMass3_[eta2]))) {
+          accept = true;
+        }  // if scoreJ > scoreI
+      }    //If high mass
+      else if (mass > lowMassCut_ && mass < highMassCut_) {
+        if (mvaScorei >= mvaScorej && ((mvaScorei > leadCutLowMass1_[eta1] && mvaScorej > subCutLowMass1_[eta2]) ||
+                                       (mvaScorei > leadCutLowMass2_[eta1] && mvaScorej > subCutLowMass2_[eta2]) ||
+                                       (mvaScorei > leadCutLowMass3_[eta1] && mvaScorej > subCutLowMass3_[eta2]))) {
+          accept = true;
+        }  //if scoreI > scoreJ
+        else if (mvaScorej > mvaScorei && ((mvaScorej > leadCutLowMass1_[eta1] && mvaScorei > subCutLowMass1_[eta2]) ||
+                                           (mvaScorej > leadCutLowMass2_[eta1] && mvaScorei > subCutLowMass2_[eta2]) ||
+                                           (mvaScorej > leadCutLowMass3_[eta1] && mvaScorei > subCutLowMass3_[eta2]))) {
+          accept = true;
+        }  //if scoreJ > scoreI
+      }    //If low mass
+    }      //j loop
+  }        //i loop
+  return accept;
+}  //Definition
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(HLTEgammaDoubleXGBoostCombFilter);

RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTClusterShapeProducer.cc

@@ -63,6 +63,7 @@ EgammaHLTClusterShapeProducer::EgammaHLTClusterShapeProducer(const edm::Paramete
   produces<reco::RecoEcalCandidateIsolationMap>("sigmaIPhiIPhi5x5NoiseCleaned");
   produces<reco::RecoEcalCandidateIsolationMap>("sMajor");
   produces<reco::RecoEcalCandidateIsolationMap>("sMinor");
+  produces<reco::RecoEcalCandidateIsolationMap>("e2x2");
 }
 
 EgammaHLTClusterShapeProducer::~EgammaHLTClusterShapeProducer() {}
@@ -108,6 +109,8 @@ void EgammaHLTClusterShapeProducer::produce(edm::StreamID sid,
   reco::RecoEcalCandidateIsolationMap clshSMajorMap(recoecalcandHandle);
   reco::RecoEcalCandidateIsolationMap clshSMinorMap(recoecalcandHandle);
 
+  reco::RecoEcalCandidateIsolationMap e2x2Map(recoecalcandHandle);
+
   for (unsigned int iRecoEcalCand = 0; iRecoEcalCand < recoecalcandHandle->size(); iRecoEcalCand++) {
     reco::RecoEcalCandidateRef recoecalcandref(recoecalcandHandle, iRecoEcalCand);
     if (recoecalcandref->superCluster()->seed()->seed().det() != DetId::Ecal) {  //HGCAL, skip for now
@@ -122,6 +125,8 @@ void EgammaHLTClusterShapeProducer::produce(edm::StreamID sid,
       clshSMajorMap.insert(recoecalcandref, 0);
       clshSMinorMap.insert(recoecalcandref, 0);
 
+      e2x2Map.insert(recoecalcandref, 0);
+
       continue;
     }
 
@@ -161,6 +166,9 @@ void EgammaHLTClusterShapeProducer::produce(edm::StreamID sid,
     float sMin = moments.sMin;
     clshSMajorMap.insert(recoecalcandref, sMaj);
     clshSMinorMap.insert(recoecalcandref, sMin);
+
+    auto const e2x2 = lazyTools.e2x2(*(recoecalcandref->superCluster()->seed()));
+    e2x2Map.insert(recoecalcandref, e2x2);
   }
 
   iEvent.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(clshMap));
@@ -175,6 +183,8 @@ void EgammaHLTClusterShapeProducer::produce(edm::StreamID sid,
 
   iEvent.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(clshSMajorMap), "sMajor");
   iEvent.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(clshSMinorMap), "sMinor");
+
+  iEvent.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(e2x2Map), "e2x2");
 }
 
 #include "FWCore/Framework/interface/MakerMacros.h"

RecoEgamma/PhotonIdentification/BuildFile.xml

@@ -22,6 +22,7 @@
 <use name="RecoLocalCalo/HcalRecAlgos"/>
 <use name="PhysicsTools/TensorFlow" />
 <use name="CommonTools/MVAUtils" />
+<use name="xgboost"/>
 <export>
   <lib name="1"/>
 </export>

RecoEgamma/PhotonIdentification/interface/PhotonXGBoostEstimator.h

@@ -0,0 +1,28 @@
+#ifndef ReciEgamma_PhotonIdentification_PhotonXGBoostEstimator_h
+#define ReciEgamma_PhotonIdentification_PhotonXGBoostEstimator_h
+
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+#include "xgboost/c_api.h"
+
+class PhotonXGBoostEstimator {
+public:
+  PhotonXGBoostEstimator(const edm::FileInPath& weightsFile, int best_ntree_limit);
+  ~PhotonXGBoostEstimator();
+
+  float computeMva(float rawEnergyIn,
+                   float r9In,
+                   float sigmaIEtaIEtaIn,
+                   float etaWidthIn,
+                   float phiWidthIn,
+                   float e2x2In,
+                   float etaIn,
+                   float hOvrEIn,
+                   float ecalPFIsoIn) const;
+
+private:
+  BoosterHandle booster_;
+  int best_ntree_limit_ = -1;
+  std::string config_;
+};
+
+#endif

RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

@@ -0,0 +1,137 @@
+#include "DataFormats/Common/interface/AssociationMap.h"
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
+#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/Framework/interface/one/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/FileInPath.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "RecoEgamma/PhotonIdentification/interface/PhotonXGBoostEstimator.h"
+
+#include <memory>
+#include <vector>
+
+class PhotonXGBoostProducer : public edm::global::EDProducer<> {
+public:
+  explicit PhotonXGBoostProducer(edm::ParameterSet const&);
+  ~PhotonXGBoostProducer() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
+
+  const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> candToken_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenR9_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenHoE_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenSigmaiEtaiEta_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenE2x2_;
+  const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenIso_;
+  const edm::FileInPath mvaFileXgbB_;
+  const edm::FileInPath mvaFileXgbE_;
+  const unsigned mvaNTreeLimitB_;
+  const unsigned mvaNTreeLimitE_;
+  const double mvaThresholdEt_;
+  std::unique_ptr<PhotonXGBoostEstimator> mvaEstimatorB_;
+  std::unique_ptr<PhotonXGBoostEstimator> mvaEstimatorE_;
+};
+
+PhotonXGBoostProducer::PhotonXGBoostProducer(edm::ParameterSet const& config)
+    : candToken_(consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("candTag"))),
+      tokenR9_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagR9"))),
+      tokenHoE_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagHoE"))),
+      tokenSigmaiEtaiEta_(
+          consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagSigmaiEtaiEta"))),
+      tokenE2x2_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagE2x2"))),
+      tokenIso_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagIso"))),
+      mvaFileXgbB_(config.getParameter<edm::FileInPath>("mvaFileXgbB")),
+      mvaFileXgbE_(config.getParameter<edm::FileInPath>("mvaFileXgbE")),
+      mvaNTreeLimitB_(config.getParameter<unsigned int>("mvaNTreeLimitB")),
+      mvaNTreeLimitE_(config.getParameter<unsigned int>("mvaNTreeLimitE")),
+      mvaThresholdEt_(config.getParameter<double>("mvaThresholdEt")) {
+  mvaEstimatorB_ = std::make_unique<PhotonXGBoostEstimator>(mvaFileXgbB_, mvaNTreeLimitB_);
+  mvaEstimatorE_ = std::make_unique<PhotonXGBoostEstimator>(mvaFileXgbE_, mvaNTreeLimitE_);
+  produces<reco::RecoEcalCandidateIsolationMap>();
+}
+
+void PhotonXGBoostProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("candTag", edm::InputTag("hltEgammaCandidatesUnseeded"));
+  desc.add<edm::InputTag>("inputTagR9", edm::InputTag("hltEgammaR9IDUnseeded", "r95x5"));
+  desc.add<edm::InputTag>("inputTagHoE", edm::InputTag("hltEgammaHoverEUnseeded"));
+  desc.add<edm::InputTag>("inputTagSigmaiEtaiEta",
+                          edm::InputTag("hltEgammaClusterShapeUnseeded", "sigmaIEtaIEta5x5NoiseCleaned"));
+  desc.add<edm::InputTag>("inputTagE2x2", edm::InputTag("hltEgammaClusterShapeUnseeded", "e2x2"));
+  desc.add<edm::InputTag>("inputTagIso", edm::InputTag("hltEgammaEcalPFClusterIsoUnseeded"));
+  desc.add<edm::FileInPath>(
+      "mvaFileXgbB", edm::FileInPath("RecoEgamma/PhotonIdentification/data/XGBoost/Photon_NTL_168_Barrel_v1.bin"));
+  desc.add<edm::FileInPath>(
+      "mvaFileXgbE", edm::FileInPath("RecoEgamma/PhotonIdentification/data/XGBoost/Photon_NTL_158_Endcap_v1.bin"));
+  desc.add<unsigned int>("mvaNTreeLimitB", 168);
+  desc.add<unsigned int>("mvaNTreeLimitE", 158);
+  desc.add<double>("mvaThresholdEt", 0);
+  descriptions.addWithDefaultLabel(desc);
+}
+
+void PhotonXGBoostProducer::produce(edm::StreamID, edm::Event& event, edm::EventSetup const& setup) const {
+  const auto& recCollection = event.getHandle(candToken_);
+
+  //get hold of r9 association map
+  const auto& r9Map = event.getHandle(tokenR9_);
+
+  //get hold of HoE association map
+  const auto& hoEMap = event.getHandle(tokenHoE_);
+
+  //get hold of isolated association map
+  const auto& sigmaiEtaiEtaMap = event.getHandle(tokenSigmaiEtaiEta_);
+
+  //get hold of e2x2 (s4) association map
+  const auto& e2x2Map = event.getHandle(tokenE2x2_);
+
+  //get hold of Ecal isolation association map
+  const auto& isoMap = event.getHandle(tokenIso_);
+
+  //output
+  reco::RecoEcalCandidateIsolationMap mvaScoreMap(recCollection);
+
+  for (size_t i = 0; i < recCollection->size(); i++) {
+    edm::Ref<reco::RecoEcalCandidateCollection> ref(recCollection, i);
+
+    float etaSC = ref->eta();
+
+    float scEnergy = ref->superCluster()->energy();
+    float r9 = (*r9Map).find(ref)->val;
+    float hoe = (*hoEMap).find(ref)->val / scEnergy;
+    float siEtaiEta = (*sigmaiEtaiEtaMap).find(ref)->val;
+    float e2x2 = (*e2x2Map).find(ref)->val;
+    float iso = (*isoMap).find(ref)->val;
+
+    float rawEnergy = ref->superCluster()->rawEnergy();
+    float etaW = ref->superCluster()->etaWidth();
+    float phiW = ref->superCluster()->phiWidth();
+
+    float scEt = scEnergy / cosh(etaSC);
+    if (scEt < 0.)
+      scEt = 0.; /* first and second order terms assume non-negative energies */
+
+    float xgbScore = -100.;
+    //compute only above threshold used for training and cand filter, else store negative value into the association map.
+    if (scEt >= mvaThresholdEt_) {
+      if (std::abs(etaSC) < 1.5)
+        xgbScore = mvaEstimatorB_->computeMva(rawEnergy, r9, siEtaiEta, etaW, phiW, e2x2, etaSC, hoe, iso);
+      else
+        xgbScore = mvaEstimatorE_->computeMva(rawEnergy, r9, siEtaiEta, etaW, phiW, e2x2, etaSC, hoe, iso);
+    }
+    mvaScoreMap.insert(ref, xgbScore);
+  }
+  event.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(mvaScoreMap));
+}
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(PhotonXGBoostProducer);