Skip to content

Commit

Permalink
Merge pull request #5 from Sam-Harper/HGCal_ShowerShapes_11_2_0_pre7_…
Browse files Browse the repository at this point in the history 
…WithExample

Adding Example config for Reco:

cmsDriver.py hltTest --python_filename egHLTPhaseIITest_cfg.py --eventcontent AODSIM --customise SLHCUpgradeSimulations/Configuration/aging.customise_aging_1000,Configuration/DataProcessing/Utils.addMonitoring --step RAW2DIGI,RECO --datatier AODSIM  --conditions auto:phase2_realistic_T15 --geometry Extended2026D49 --nThreads 4 --filein /store/mc/Phase2HLTTDRSummer20ReRECOMiniAOD/DoubleElectron_FlatPt-1To100/GEN-SIM-DIGI-RAW-MINIAOD/PU200_111X_mcRun4_realistic_T15_v1-v2/280000/1B09CC5D-FAA3-5442-BE7B-96302835E860.root --era Phase2C9 --no_exec --mc --customise_commands "process.load('RecoEgamma.EgammaHLTProducers.HLTEgPhaseIITestSequence_cff');process.HLT_EgPhaseIITest_v1 = cms.Path(process.HLTEgPhaseIITestSequence);process.schedule.append(process.HLT_EgPhaseIITest_v1);process.AODSIMoutput.outputCommands.extend(['keep recoEcalCandidates*_*_*_*','keep recoRecoEcalCandidatesToValuefloatAssociation_*_*_*'])" -n 12 --processName=RECO2
  • Loading branch information
@SohamBhattacharya
SohamBhattacharya committed Jan 5, 2021
2 parents 58c99c2 + 833d51f commit 994d7f9
Show file tree
Hide file tree
Showing 5 changed files with 2,882 additions and 34 deletions.
15 changes: 8 additions & 7 deletions RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc
View file
Expand Up @@ -54,6 +54,7 @@ class EgammaHLTHGCalIDVarProducer : public edm::stream::EDProducer<> {
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)
Expand All @@ -62,7 +63,9 @@ EgammaHLTHGCalIDVarProducer::EgammaHLTHGCalIDVarProducer(const edm::ParameterSet
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"))) {
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"));
Expand Down Expand Up @@ -97,9 +100,7 @@ void EgammaHLTHGCalIDVarProducer::produce(edm::Event& iEvent, const edm::EventSe
const auto& hgcalRecHits = iEvent.get(hgcalRecHitToken_);
const auto& layerClusters = iEvent.get(layerClusterToken_);

auto theconsumes = consumesCollector();
HGCalShowerShapeHelper ssCalc(theconsumes);
ssCalc.initPerEvent(iSetup, hgcalRecHits);
ssCalc_.initPerEvent(iSetup, hgcalRecHits);

auto rVarMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
auto hForHoverEMap = std::make_unique<reco::RecoEcalCandidateIsolationMap>(recoEcalCandHandle);
Expand All @@ -109,13 +110,13 @@ void EgammaHLTHGCalIDVarProducer::produce(edm::Event& iEvent, const edm::EventSe

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()));
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_);
auto pcaWidths = ssCalc_.getPCAWidths(rCylinder_);
for (auto& pcaMap : pcaAssocMaps_) {
pcaMap.insert(candRef, pcaWidths);
}
Expand Down
190 changes: 190 additions & 0 deletions RecoEgamma/EgammaHLTProducers/plugins/HLTHGCalLayerClusterIsolationProducer.cc
View file
@@ -0,0 +1,190 @@
#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);

0 comments on commit 994d7f9

Please sign in to comment.