Skip to content

Commit

Permalink
Merge pull request #32520 from SohamBhattacharya/HGCal_ShowerShapes_1…
Browse files Browse the repository at this point in the history
…1_2_0_pre7

HGCal ID variables
  • Loading branch information
cmsbuild authored Jan 11, 2021
2 parents b3e3916 + de4af6a commit e06cb88
Show file tree
Hide file tree
Showing 8 changed files with 3,506 additions and 0 deletions.
1 change: 1 addition & 0 deletions RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
<use name="Geometry/CaloGeometry"/>
<use name="RecoEcal/EgammaCoreTools"/>
<use name="RecoEgamma/EgammaIsolationAlgos"/>
<use name="RecoEgamma/EgammaTools"/>
<use name="MagneticField/Engine"/>
<use name="MagneticField/Records"/>
<use name="DataFormats/DetId"/>
Expand Down
130 changes: 130 additions & 0 deletions RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h"
#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"

class EgammaHLTHGCalIDVarProducer : public edm::stream::EDProducer<> {
public:
explicit EgammaHLTHGCalIDVarProducer(const edm::ParameterSet&);
~EgammaHLTHGCalIDVarProducer() override;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
void produce(edm::Event&, const edm::EventSetup&) override;

class PCAAssocMap {
public:
PCAAssocMap(double HGCalShowerShapeHelper::ShowerWidths::*var, const std::string& name) : var_(var), name_(name) {}

void initMap(const edm::Handle<reco::RecoEcalCandidateCollection>& candHandle) {
assocMap_ = std::make_unique<reco::RecoEcalCandidateIsolationMap>(candHandle);
}

void insert(reco::RecoEcalCandidateRef& ref, const HGCalShowerShapeHelper::ShowerWidths& showerWidths) {
assocMap_->insert(ref, showerWidths.*var_);
}

std::unique_ptr<reco::RecoEcalCandidateIsolationMap> releaseMap() { return std::move(assocMap_); }
const std::string& name() const { return name_; }

private:
double HGCalShowerShapeHelper::ShowerWidths::*var_;
std::string name_;
std::unique_ptr<reco::RecoEcalCandidateIsolationMap> assocMap_;
};

private:
// ----------member data ---------------------------
float rCylinder_;
float hOverECone_;
std::vector<PCAAssocMap> pcaAssocMaps_;
const edm::EDGetTokenT<reco::RecoEcalCandidateCollection> recoEcalCandidateToken_;
const edm::EDGetTokenT<reco::PFRecHitCollection> hgcalRecHitToken_;
const edm::EDGetTokenT<reco::CaloClusterCollection> layerClusterToken_;
HGCalShowerShapeHelper ssCalc_;
};

EgammaHLTHGCalIDVarProducer::EgammaHLTHGCalIDVarProducer(const edm::ParameterSet& config)
: rCylinder_(config.getParameter<double>("rCylinder")),
hOverECone_(config.getParameter<double>("hOverECone")),
recoEcalCandidateToken_(
consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
hgcalRecHitToken_(consumes<reco::PFRecHitCollection>(config.getParameter<edm::InputTag>("hgcalRecHits"))),
layerClusterToken_(consumes<reco::CaloClusterCollection>(config.getParameter<edm::InputTag>("layerClusters"))),
ssCalc_(consumesCollector()) {
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xx, "sigma2xx"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yy, "sigma2yy"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zz, "sigma2zz"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xy, "sigma2xy"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yz, "sigma2yz"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zx, "sigma2zx"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2uu, "sigma2uu"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2vv, "sigma2vv"));
pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2ww, "sigma2ww"));

produces<reco::RecoEcalCandidateIsolationMap>("rVar");
produces<reco::RecoEcalCandidateIsolationMap>("hForHOverE");
for (auto& var : pcaAssocMaps_) {
produces<reco::RecoEcalCandidateIsolationMap>(var.name());
}
}

EgammaHLTHGCalIDVarProducer::~EgammaHLTHGCalIDVarProducer() {}

void EgammaHLTHGCalIDVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("recoEcalCandidateProducer", edm::InputTag("hltL1SeededRecoEcalCandidate"));
desc.add<edm::InputTag>("hgcalRecHits", edm::InputTag("hgcalRecHits"));
desc.add<edm::InputTag>("layerClusters", edm::InputTag("layerClusters"));
desc.add<double>("rCylinder", 2.8);
desc.add<double>("hOverECone", 0.15);
descriptions.add(("hltEgammaHLTHGCalIDVarProducer"), desc);
}

void EgammaHLTHGCalIDVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateToken_);
const auto& hgcalRecHits = iEvent.get(hgcalRecHitToken_);
const auto& layerClusters = iEvent.get(layerClusterToken_);

ssCalc_.initPerEvent(iSetup, hgcalRecHits);

auto rVarMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
auto hForHoverEMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
for (auto& pcaMap : pcaAssocMaps_) {
pcaMap.initMap(recoEcalCandHandle);
}

for (size_t candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) {
reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr);
ssCalc_.initPerObject(candRef->superCluster()->hitsAndFractions());
rVarMap->insert(candRef, ssCalc_.getRvar(rCylinder_, candRef->superCluster()->energy()));

float hForHoverE = HGCalClusterTools::hadEnergyInCone(
candRef->superCluster()->eta(), candRef->superCluster()->phi(), layerClusters, 0., hOverECone_, 0., 0.);
hForHoverEMap->insert(candRef, hForHoverE);
auto pcaWidths = ssCalc_.getPCAWidths(rCylinder_);
for (auto& pcaMap : pcaAssocMaps_) {
pcaMap.insert(candRef, pcaWidths);
}
}
iEvent.put(std::move(rVarMap), "rVar");
iEvent.put(std::move(hForHoverEMap), "hForHOverE");
for (auto& pcaMap : pcaAssocMaps_) {
iEvent.put(pcaMap.releaseMap(), pcaMap.name());
}
}

DEFINE_FWK_MODULE(EgammaHLTHGCalIDVarProducer);
Original file line number Diff line number Diff line change
@@ -0,0 +1,197 @@
#include <iostream>
#include <vector>
#include <memory>

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"

#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateIsolation.h"

#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"

#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"

template <typename T1>
class HLTHGCalLayerClusterIsolationProducer : public edm::stream::EDProducer<> {
typedef std::vector<T1> T1Collection;
typedef edm::Ref<T1Collection> T1Ref;
typedef edm::AssociationMap<edm::OneToValue<std::vector<T1>, float>> T1IsolationMap;

public:
explicit HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet&);
~HLTHGCalLayerClusterIsolationProducer() override = default;

void produce(edm::Event&, const edm::EventSetup&) override;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
edm::EDGetTokenT<T1Collection> recoCandidateProducer_;
const edm::EDGetTokenT<reco::CaloClusterCollection> layerClusterProducer_;
const edm::EDGetTokenT<double> rhoProducer_;

const double drMax_;
const double drVetoEM_;
const double drVetoHad_;
const double minEnergyEM_;
const double minEnergyHad_;
const double minEtEM_;
const double minEtHad_;
const bool useEt_;
const bool doRhoCorrection_;
const double rhoMax_;
const double rhoScale_;
const std::vector<double> effectiveAreas_;
};

template <typename T1>
HLTHGCalLayerClusterIsolationProducer<T1>::HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet& config)
: layerClusterProducer_(
consumes<reco::CaloClusterCollection>(config.getParameter<edm::InputTag>("layerClusterProducer"))),
rhoProducer_(consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))),
drMax_(config.getParameter<double>("drMax")),
drVetoEM_(config.getParameter<double>("drVetoEM")),
drVetoHad_(config.getParameter<double>("drVetoHad")),
minEnergyEM_(config.getParameter<double>("minEnergyEM")),
minEnergyHad_(config.getParameter<double>("minEnergyHad")),
minEtEM_(config.getParameter<double>("minEtEM")),
minEtHad_(config.getParameter<double>("minEtHad")),
useEt_(config.getParameter<bool>("useEt")),
doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
rhoMax_(config.getParameter<double>("rhoMax")),
rhoScale_(config.getParameter<double>("rhoScale")),
effectiveAreas_(config.getParameter<std::vector<double>>("effectiveAreas")) {
if (doRhoCorrection_) {
if (effectiveAreas_.size() != 2)
throw cms::Exception("IncompatibleVects")
<< "effectiveAreas should have two elements for em and had components. \n";
}

std::string recoCandidateProducerName = "recoCandidateProducer";
if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
recoCandidateProducerName = "recoEcalCandidateProducer";

recoCandidateProducer_ = consumes<T1Collection>(config.getParameter<edm::InputTag>(recoCandidateProducerName));
produces<T1IsolationMap>();
produces<T1IsolationMap>("em");
produces<T1IsolationMap>("had");
}

template <typename T1>
void HLTHGCalLayerClusterIsolationProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
std::string recoCandidateProducerName = "recoCandidateProducer";
if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
recoCandidateProducerName = "recoEcalCandidateProducer";

edm::ParameterSetDescription desc;
desc.add<edm::InputTag>(recoCandidateProducerName, edm::InputTag("hltL1SeededRecoEcalCandidatePF"));
desc.add<edm::InputTag>("layerClusterProducer", edm::InputTag("hltParticleFlowClusterECAL"));
desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
desc.add<bool>("doRhoCorrection", false);
desc.add<bool>("useEt", false);
desc.add<double>("rhoMax", 9.9999999E7);
desc.add<double>("rhoScale", 1.0);
desc.add<double>("drMax", 0.3);
desc.add<double>("drVetoEM", 0.0);
desc.add<double>("drVetoHad", 0.0);
desc.add<double>("minEnergyEM", 0.0);
desc.add<double>("minEnergyHad", 0.0);
desc.add<double>("minEtEM", 0.0);
desc.add<double>("minEtHad", 0.0);
desc.add<std::vector<double>>("effectiveAreas", {0.0, 0.0}); // for em and had components
descriptions.add(defaultModuleLabel<HLTHGCalLayerClusterIsolationProducer<T1>>(), desc);
}

template <typename T1>
void HLTHGCalLayerClusterIsolationProducer<T1>::produce(edm::Event& iEvent, const edm::EventSetup&) {
edm::Handle<double> rhoHandle;
double rho = 0.0;
if (doRhoCorrection_) {
iEvent.getByToken(rhoProducer_, rhoHandle);
rho = *(rhoHandle.product());
}

rho = std::min(rho, rhoMax_);
rho = rho * rhoScale_;

edm::Handle<T1Collection> recoCandHandle;
edm::Handle<reco::CaloClusterCollection> clusterHandle;

iEvent.getByToken(recoCandidateProducer_, recoCandHandle);
iEvent.getByToken(layerClusterProducer_, clusterHandle);

const std::vector<reco::CaloCluster> layerClusters = *(clusterHandle.product());

T1IsolationMap recoCandMap(recoCandHandle);
T1IsolationMap recoCandMapEm(recoCandHandle);
T1IsolationMap recoCandMapHad(recoCandHandle);

for (unsigned int iReco = 0; iReco < recoCandHandle->size(); iReco++) {
T1Ref candRef(recoCandHandle, iReco);

float sumEm =
HGCalClusterTools::emEnergyInCone(candRef->eta(),
candRef->phi(),
layerClusters,
drVetoEM_,
drMax_,
minEtEM_,
minEnergyEM_,
useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);

float sumHad =
HGCalClusterTools::hadEnergyInCone(candRef->eta(),
candRef->phi(),
layerClusters,
drVetoHad_,
drMax_,
minEtHad_,
minEnergyHad_,
useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);

if (doRhoCorrection_) {
sumEm = sumEm - rho * effectiveAreas_.at(0);
sumHad = sumHad - rho * effectiveAreas_.at(1);
}

float sum = sumEm + sumHad;

recoCandMap.insert(candRef, sum);
recoCandMapEm.insert(candRef, sumEm);
recoCandMapHad.insert(candRef, sumHad);
}

iEvent.put(std::make_unique<T1IsolationMap>(recoCandMap));
iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapEm), "em");
iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapHad), "had");
}

typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate> EgammaHLTHGCalLayerClusterIsolationProducer;
typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoChargedCandidate> MuonHLTHGCalLayerClusterIsolationProducer;

DEFINE_FWK_MODULE(EgammaHLTHGCalLayerClusterIsolationProducer);
DEFINE_FWK_MODULE(MuonHLTHGCalLayerClusterIsolationProducer);
Loading

0 comments on commit e06cb88

Please sign in to comment.