Skip to content

Commit

Permalink
* e2x2 produced by the EGammaClusterShapeProducer
Browse files Browse the repository at this point in the history
* Photon XGBoost MVA estimator and producer.
Implemented using native XGBoost C API for photon MVA score
Currently diphoton score filter is not yet submitted.

Model files will be added to the cms-data

XGBoost Double photon combination filter

fix module name

code-format

fix comments from maintainers

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc

Co-authored-by: Marco Musich <[email protected]>

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>

tidy up and sort headers

code-format

Update RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc

Co-authored-by: Marco Musich <[email protected]>
  • Loading branch information
smorovic and mmusich committed Mar 21, 2024
1 parent c374ae0 commit b2bdbef
Show file tree
Hide file tree
Showing 6 changed files with 393 additions and 0 deletions.
157 changes: 157 additions & 0 deletions HLTrigger/Egamma/plugins/HLTEgammaDoubleXGBoostCombFilter.cc
Original file line number Diff line number Diff line change
@@ -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);
Original file line number Diff line number Diff line change
Expand Up @@ -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() {}
Expand Down Expand Up @@ -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
Expand All @@ -122,6 +125,8 @@ void EgammaHLTClusterShapeProducer::produce(edm::StreamID sid,
clshSMajorMap.insert(recoecalcandref, 0);
clshSMinorMap.insert(recoecalcandref, 0);

e2x2Map.insert(recoecalcandref, 0);

continue;
}

Expand Down Expand Up @@ -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));
Expand All @@ -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"
Expand Down
1 change: 1 addition & 0 deletions RecoEgamma/PhotonIdentification/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
<use name="RecoLocalCalo/HcalRecAlgos"/>
<use name="PhysicsTools/TensorFlow" />
<use name="CommonTools/MVAUtils" />
<use name="xgboost"/>
<export>
<lib name="1"/>
</export>
28 changes: 28 additions & 0 deletions RecoEgamma/PhotonIdentification/interface/PhotonXGBoostEstimator.h
Original file line number Diff line number Diff line change
@@ -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
137 changes: 137 additions & 0 deletions RecoEgamma/PhotonIdentification/plugins/PhotonXGBoostProducer.cc
Original file line number Diff line number Diff line change
@@ -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/xgb_photonmva_barrel_v1.bin"));
desc.add<edm::FileInPath>("mvaFileXgbE",
edm::FileInPath("RecoEgamma/PhotonIdentification/data/xgb_photonmva_endcap_v1.bin"));
desc.add<unsigned int>("mvaNTreeLimitB", 55);
desc.add<unsigned int>("mvaNTreeLimitE", 48);
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);
Loading

0 comments on commit b2bdbef

Please sign in to comment.