PhotonIdentification/plugins/PhotonXGBoostProducer.cc

0001 #include "DataFormats/Common/interface/AssociationMap.h"
0002 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0003 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0004 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0005 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007 #include "FWCore/Framework/interface/one/EDProducer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0011 #include "FWCore/ParameterSet/interface/FileInPath.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "RecoEgamma/PhotonIdentification/interface/PhotonXGBoostEstimator.h"
0016
0017 #include <memory>
0018 #include <vector>
0019
0020 class PhotonXGBoostProducer : public edm::global::EDProducer<> {
0021 public:
0022   explicit PhotonXGBoostProducer(edm::ParameterSet const&);
0023   ~PhotonXGBoostProducer() = default;
0024
0025   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0026
0027 private:
0028   void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
0029
0030   const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> candToken_;
0031   const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenR9_;
0032   const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenHoE_;
0033   const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenSigmaiEtaiEta_;
0034   const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenE2x2_;
0035   const edm::EDGetTokenT<reco::RecoEcalCandidateIsolationMap> tokenIso_;
0036   const edm::FileInPath mvaFileXgbB_;
0037   const edm::FileInPath mvaFileXgbE_;
0038   const unsigned mvaNTreeLimitB_;
0039   const unsigned mvaNTreeLimitE_;
0040   const double mvaThresholdEt_;
0041   const std::unique_ptr<const PhotonXGBoostEstimator> mvaEstimatorB_;
0042   const std::unique_ptr<const PhotonXGBoostEstimator> mvaEstimatorE_;
0043 };
0044
0045 PhotonXGBoostProducer::PhotonXGBoostProducer(edm::ParameterSet const& config)
0046     : candToken_(consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("candTag"))),
0047       tokenR9_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagR9"))),
0048       tokenHoE_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagHoE"))),
0049       tokenSigmaiEtaiEta_(
0050           consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagSigmaiEtaiEta"))),
0051       tokenE2x2_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagE2x2"))),
0052       tokenIso_(consumes<reco::RecoEcalCandidateIsolationMap>(config.getParameter<edm::InputTag>("inputTagIso"))),
0053       mvaFileXgbB_(config.getParameter<edm::FileInPath>("mvaFileXgbB")),
0054       mvaFileXgbE_(config.getParameter<edm::FileInPath>("mvaFileXgbE")),
0055       mvaNTreeLimitB_(config.getParameter<unsigned int>("mvaNTreeLimitB")),
0056       mvaNTreeLimitE_(config.getParameter<unsigned int>("mvaNTreeLimitE")),
0057       mvaThresholdEt_(config.getParameter<double>("mvaThresholdEt")),
0058       mvaEstimatorB_{std::make_unique<const PhotonXGBoostEstimator>(mvaFileXgbB_, mvaNTreeLimitB_)},
0059       mvaEstimatorE_{std::make_unique<const PhotonXGBoostEstimator>(mvaFileXgbE_, mvaNTreeLimitE_)} {
0060   produces<reco::RecoEcalCandidateIsolationMap>();
0061 }
0062
0063 void PhotonXGBoostProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0064   edm::ParameterSetDescription desc;
0065   desc.add<edm::InputTag>("candTag", edm::InputTag("hltEgammaCandidatesUnseeded"));
0066   desc.add<edm::InputTag>("inputTagR9", edm::InputTag("hltEgammaR9IDUnseeded", "r95x5"));
0067   desc.add<edm::InputTag>("inputTagHoE", edm::InputTag("hltEgammaHoverEUnseeded"));
0068   desc.add<edm::InputTag>("inputTagSigmaiEtaiEta",
0069                           edm::InputTag("hltEgammaClusterShapeUnseeded", "sigmaIEtaIEta5x5NoiseCleaned"));
0070   desc.add<edm::InputTag>("inputTagE2x2", edm::InputTag("hltEgammaClusterShapeUnseeded", "e2x2"));
0071   desc.add<edm::InputTag>("inputTagIso", edm::InputTag("hltEgammaEcalPFClusterIsoUnseeded"));
0072   desc.add<edm::FileInPath>(
0073       "mvaFileXgbB", edm::FileInPath("RecoEgamma/PhotonIdentification/data/XGBoost/Photon_NTL_168_Barrel_v1.bin"));
0074   desc.add<edm::FileInPath>(
0075       "mvaFileXgbE", edm::FileInPath("RecoEgamma/PhotonIdentification/data/XGBoost/Photon_NTL_158_Endcap_v1.bin"));
0076   desc.add<unsigned int>("mvaNTreeLimitB", 168);
0077   desc.add<unsigned int>("mvaNTreeLimitE", 158);
0078   desc.add<double>("mvaThresholdEt", 0);
0079   descriptions.addWithDefaultLabel(desc);
0080 }
0081
0082 void PhotonXGBoostProducer::produce(edm::StreamID, edm::Event& event, edm::EventSetup const& setup) const {
0083   const auto& recCollection = event.getHandle(candToken_);
0084
0085   //get hold of r9 association map
0086   const auto& r9Map = event.getHandle(tokenR9_);
0087
0088   //get hold of HoE association map
0089   const auto& hoEMap = event.getHandle(tokenHoE_);
0090
0091   //get hold of isolated association map
0092   const auto& sigmaiEtaiEtaMap = event.getHandle(tokenSigmaiEtaiEta_);
0093
0094   //get hold of e2x2 (s4) association map
0095   const auto& e2x2Map = event.getHandle(tokenE2x2_);
0096
0097   //get hold of Ecal isolation association map
0098   const auto& isoMap = event.getHandle(tokenIso_);
0099
0100   //output
0101   reco::RecoEcalCandidateIsolationMap mvaScoreMap(recCollection);
0102
0103   for (size_t i = 0; i < recCollection->size(); i++) {
0104     edm::Ref<reco::RecoEcalCandidateCollection> ref(recCollection, i);
0105
0106     float etaSC = ref->eta();
0107
0108     float scEnergy = ref->superCluster()->energy();
0109     float r9 = (*r9Map).find(ref)->val;
0110     float hoe = (*hoEMap).find(ref)->val / scEnergy;
0111     float siEtaiEta = (*sigmaiEtaiEtaMap).find(ref)->val;
0112     float e2x2 = (*e2x2Map).find(ref)->val;
0113     float s4 = e2x2 / scEnergy;
0114     float iso = (*isoMap).find(ref)->val;
0115
0116     float rawEnergy = ref->superCluster()->rawEnergy();
0117     float etaW = ref->superCluster()->etaWidth();
0118     float phiW = ref->superCluster()->phiWidth();
0119
0120     float scEt = scEnergy / cosh(etaSC);
0121     if (scEt < 0.)
0122       scEt = 0.; /* first and second order terms assume non-negative energies */
0123
0124     float xgbScore = -100.;
0125     //compute only above threshold used for training and cand filter, else store negative value into the association map.
0126     if (scEt >= mvaThresholdEt_) {
0127       if (std::abs(etaSC) < 1.5)
0128         xgbScore = mvaEstimatorB_->computeMva(rawEnergy, r9, siEtaiEta, etaW, phiW, s4, etaSC, hoe, iso);
0129       else
0130         xgbScore = mvaEstimatorE_->computeMva(rawEnergy, r9, siEtaiEta, etaW, phiW, s4, etaSC, hoe, iso);
0131     }
0132     mvaScoreMap.insert(ref, xgbScore);
0133   }
0134   event.put(std::make_unique<reco::RecoEcalCandidateIsolationMap>(mvaScoreMap));
0135 }
0136
0137 #include "FWCore/Framework/interface/MakerMacros.h"
0138 DEFINE_FWK_MODULE(PhotonXGBoostProducer);