-
Notifications
You must be signed in to change notification settings - Fork 4.2k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
HGCal ID variables #32519
HGCal ID variables #32519
Changes from 10 commits
41f3043
89a2e50
3f9799d
2f8be5a
9da375d
9afea05
bf57b1d
5fdf813
e3e931f
2626651
1426f18
58c99c2
da3597c
bc80b81
02d7790
833d51f
994d7f9
5612e4d
de4af6a
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
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_; | ||
}; | ||
|
||
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"))) { | ||
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_); | ||
|
||
auto theconsumes = consumesCollector(); | ||
HGCalShowerShapeHelper ssCalc(theconsumes); | ||
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,44 @@ | ||
#ifndef RecoEgamma_EgammaTools_HGCalClusterTools_h | ||
#define RecoEgamma_EgammaTools_HGCalClusterTools_h | ||
|
||
#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" | ||
#include <vector> | ||
|
||
class HGCalClusterTools { | ||
public: | ||
enum class EType { ET, ENERGY }; | ||
|
||
static float energyInCone(const float eta, | ||
const float phi, | ||
const std::vector<reco::CaloCluster>& layerClusters, | ||
const float minDR, | ||
const float maxDR, | ||
const float minEt, | ||
const float minEnergy, | ||
const std::vector<DetId::Detector>& subDets, | ||
const HGCalClusterTools::EType& eType = EType::ENERGY); | ||
|
||
static float hadEnergyInCone(const float eta, | ||
const float phi, | ||
const std::vector<reco::CaloCluster>& layerClusters, | ||
const float minDR, | ||
const float maxDR, | ||
const float minEt, | ||
const float minEnergy, | ||
const HGCalClusterTools::EType& eType = EType::ENERGY) { | ||
return energyInCone( | ||
eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalHSi, DetId::HGCalHSc}, eType); | ||
} | ||
static float emEnergyInCone(const float eta, | ||
const float phi, | ||
const std::vector<reco::CaloCluster>& layerClusters, | ||
const float minDR, | ||
const float maxDR, | ||
const float minEt, | ||
const float minEnergy, | ||
const HGCalClusterTools::EType& eType = EType::ENERGY) { | ||
return energyInCone(eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalEE}, eType); | ||
} | ||
}; | ||
|
||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,132 @@ | ||
#ifndef RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h | ||
#define RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h | ||
|
||
// system include files | ||
#include <algorithm> | ||
#include <cstdlib> | ||
#include <iostream> | ||
#include <map> | ||
#include <memory> | ||
#include <string> | ||
#include <utility> | ||
#include <vector> | ||
|
||
// external include files | ||
#include <CLHEP/Vector/LorentzVector.h> | ||
#include <Math/Point3D.h> | ||
#include <Math/Point3Dfwd.h> | ||
#include <TMatrixD.h> | ||
#include <TVectorD.h> | ||
|
||
// CMSSW include files | ||
#include "DataFormats/CaloRecHit/interface/CaloCluster.h" | ||
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" | ||
#include "DataFormats/FWLite/interface/ESHandle.h" | ||
#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h" | ||
#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" | ||
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" | ||
#include "DataFormats/HGCRecHit/interface/HGCRecHit.h" | ||
#include "DataFormats/Math/interface/LorentzVector.h" | ||
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" | ||
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h" | ||
#include "DataFormats/TrackReco/interface/Track.h" | ||
#include "DataFormats/TrackReco/interface/TrackFwd.h" | ||
#include "FWCore/Framework/interface/Event.h" | ||
#include "FWCore/Framework/interface/ESHandle.h" | ||
#include "FWCore/Framework/interface/Frameworkfwd.h" | ||
#include "FWCore/Framework/interface/MakerMacros.h" | ||
#include "FWCore/Framework/interface/stream/EDProducer.h" | ||
#include "FWCore/ParameterSet/interface/ParameterSet.h" | ||
#include "FWCore/Utilities/interface/StreamID.h" | ||
#include "Geometry/CaloTopology/interface/HGCalTopology.h" | ||
#include "Geometry/Records/interface/IdealGeometryRecord.h" | ||
#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" | ||
#include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h" | ||
|
||
class HGCalShowerShapeHelper { | ||
private: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. public, protected, private must be declared in that order (4.25) (I missed pointing this out earlier, thanks Slava for pointing out) |
||
// Good to filter/compute/store this stuff beforehand as they are common to the shower shape variables. | ||
// No point in filtering, computing layer-wise centroids, etc. for each variable again and again. | ||
// Once intitialized, one can the calculate different variables one after another for a given object. | ||
// If a different set of preselections (E, ET, etc.) is required for a given object, then reinitialize using initPerObject(...). | ||
|
||
// In principle should consider the HGCalHSi and HGCalHSc hits (leakage) also. | ||
// Can have subdetector dependent thresholds and layer selection. | ||
// To be implemented. | ||
|
||
double minHitE_; | ||
double minHitET_; | ||
double minHitET2_; | ||
int minLayer_; | ||
int maxLayer_; | ||
int nLayer_; | ||
DetId::Detector subDet_; | ||
|
||
edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_; | ||
hgcal::RecHitTools recHitTools_; | ||
|
||
std::unordered_map<uint32_t, const reco::PFRecHit *> pfRecHitPtrMap_; | ||
std::vector<std::pair<DetId, float> > hitsAndFracs_; | ||
std::vector<double> hitEnergies_; | ||
std::vector<double> hitEnergiesWithFracs_; | ||
|
||
ROOT::Math::XYZVector centroid_; | ||
std::vector<double> layerEnergies_; | ||
std::vector<ROOT::Math::XYZVector> layerCentroids_; | ||
|
||
void setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits); | ||
|
||
void setFilteredHitsAndFractions(const std::vector<std::pair<DetId, float> > &hitsAndFracs); | ||
|
||
public: | ||
static const double kLDWaferCellSize_; | ||
static const double kHDWaferCellSize_; | ||
|
||
void setLayerWiseInfo(); | ||
|
||
struct ShowerWidths { | ||
double sigma2xx; | ||
double sigma2yy; | ||
double sigma2zz; | ||
|
||
double sigma2xy; | ||
double sigma2yz; | ||
double sigma2zx; | ||
|
||
double sigma2uu; | ||
double sigma2vv; | ||
double sigma2ww; | ||
|
||
ShowerWidths() | ||
: sigma2xx(0.0), | ||
sigma2yy(0.0), | ||
sigma2zz(0.0), | ||
sigma2xy(0.0), | ||
sigma2yz(0.0), | ||
sigma2zx(0.0), | ||
sigma2uu(0.0), | ||
sigma2vv(0.0), | ||
sigma2ww(0.0) {} | ||
}; | ||
|
||
HGCalShowerShapeHelper(edm::ConsumesCollector &sumes); | ||
|
||
void initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &recHits); | ||
|
||
void initPerObject(const std::vector<std::pair<DetId, float> > &hitsAndFracs, | ||
double minHitE = 0, | ||
double minHitET = 0, | ||
int minLayer = 1, | ||
int maxLayer = -1, | ||
DetId::Detector subDet = DetId::HGCalEE); | ||
|
||
const double getCellSize(DetId detId); | ||
|
||
// Compute Rvar in a cylinder around the layer centroids | ||
const double getRvar(double cylinderR, double energyNorm, bool useFractions = true, bool useCellSize = true); | ||
|
||
// Compute PCA widths around the layer centroids | ||
const ShowerWidths getPCAWidths(double cylinderR, bool useFractions = false); | ||
}; | ||
|
||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,56 @@ | ||
#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h" | ||
|
||
#include "DataFormats/CaloRecHit/interface/CaloCluster.h" | ||
#include "DataFormats/Math/interface/deltaR.h" | ||
#include "DataFormats/Math/interface/deltaPhi.h" | ||
|
||
float HGCalClusterTools::energyInCone(const float eta, | ||
const float phi, | ||
const std::vector<reco::CaloCluster>& layerClusters, | ||
const float minDR, | ||
const float maxDR, | ||
const float minEt, | ||
const float minEnergy, | ||
const std::vector<DetId::Detector>& subDets, | ||
const HGCalClusterTools::EType& eType) { | ||
float hadValue = 0.; | ||
|
||
const float minDR2 = minDR * minDR; | ||
const float maxDR2 = maxDR * maxDR; | ||
|
||
for (auto& clus : layerClusters) { | ||
if (clus.energy() < minEnergy) { | ||
continue; | ||
} | ||
|
||
if (std::find(subDets.begin(), subDets.end(), clus.seed().det()) == subDets.end()) { | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. what's There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The size is small (1-3) as |
||
continue; | ||
} | ||
|
||
float clusEt = clus.energy() * std::sin(clus.position().theta()); | ||
if (clusEt < minEt) { | ||
continue; | ||
} | ||
|
||
//this is a prefilter on the clusters before we calculuate | ||
//the expensive eta() of the cluster | ||
float dPhi = reco::deltaPhi(phi, clus.phi()); | ||
if (dPhi > maxDR) { | ||
continue; | ||
} | ||
|
||
float dR2 = reco::deltaR2(eta, phi, clus.eta(), clus.phi()); | ||
if (dR2 < minDR2 || dR2 > maxDR2) { | ||
continue; | ||
} | ||
switch (eType) { | ||
case EType::ET: | ||
hadValue += clusEt; | ||
break; | ||
case EType::ENERGY: | ||
hadValue += clus.energy(); | ||
break; | ||
} | ||
} | ||
return hadValue; | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The way this class is used below should be fine - but the interfaces it presents exposes makes it a bit error prone.
Returning a
std::unique_ptr<...>&
has (at least ?) two problems:unique_ptr&
is used; if it gets destroyed before the usage, one gets a dangling pointerIf the intended behaviour is that of shared ownership, it would be better to hold and return a
std::sharet_ptr
.If the intended behaviour is to transfer ownership in a well-defined place (which I think is the case here), it would be better if
initMap()
returnsvoid
(the return value is not used below)operator()
is renames to something likegetMap()
and returnsstd::unique_ptr<...>
by value, not by referenceThere was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
you're completely correct of course, will adjust