Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Adding HGCAL SuperCluster Regressions #32901

Merged
merged 6 commits into from
Feb 24, 2021
Merged

Adding HGCAL SuperCluster Regressions #32901

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
811d7e4
initial support for HGCAL in regression
Sam-Harper Jan 19, 2021
a92de71
setting phase2 GT
Sam-Harper Feb 12, 2021
42dac06
fixing phase2 reco config
Sam-Harper Feb 12, 2021
affd483
bug fixes for HLT reg
Sam-Harper Feb 12, 2021
36e6d29
code review comments
Sam-Harper Feb 19, 2021
c8f2356
removing dictionary and setting final GT
Sam-Harper Feb 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion Configuration/AlCa/python/autoCond.py
View file
Open in desktop
Expand Up @@ -83,7 +83,7 @@
# GlobalTag for MC production with realistic conditions for Phase1 2024
'phase1_2024_realistic' : '113X_mcRun3_2024_realistic_v4', # GT containing realistic conditions for Phase1 2024
# GlobalTag for MC production with realistic conditions for Phase2
'phase2_realistic' : '113X_mcRun4_realistic_v3'
'phase2_realistic' : '113X_mcRun4_realistic_v4'
}

aliases = {
Expand Down
2 changes: 2 additions & 0 deletions DataFormats/Common/src/classes_def.xml
View file
Open in desktop
Expand Up @@ -151,12 +151,14 @@
<class name="edm::ValueMap<float>" />
<class name="edm::ValueMap<double>" />
<class name="edm::ValueMap<std::pair<float, float>>" />
<class name="edm::ValueMap<std::vector<float>>" />
<class name="edm::Wrapper<edm::ValueMap<int> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<bool> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<unsigned int> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<float> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<double> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<std::pair<float, float>> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<std::vector<float>> >" splitLevel="0"/>
<class name="std::vector<edm::EventAuxiliary>"/>
<class name="edm::Wrapper<std::vector<edm::EventAuxiliary> >"/>
<class name="edm::ErrorSummaryEntry" ClassVersion="10">
Expand Down
3 changes: 3 additions & 0 deletions RecoEcal/EgammaClusterAlgos/BuildFile.xml
View file
Open in desktop
Expand Up @@ -3,12 +3,15 @@
<use name="DataFormats/EcalRecHit"/>
<use name="DataFormats/EgammaReco"/>
<use name="DataFormats/Math"/>
<use name="DataFormats/ParticleFlowReco"/>
<use name="DataFormats/VertexReco"/>
<use name="CondFormats/ESObjects"/>
<use name="CondFormats/GBRForest"/>
<use name="CondFormats/DataRecord"/>
<use name="RecoEcal/EgammaCoreTools"/>
<use name="Geometry/CaloGeometry"/>
<use name="Geometry/CaloTopology"/>
<use name="RecoEgamma/EgammaTools"/>
<use name="RecoParticleFlow/PFClusterTools"/>
<use name="clhep"/>
<export>
Expand Down
184 changes: 140 additions & 44 deletions RecoEcal/EgammaClusterAlgos/interface/SCEnergyCorrectorSemiParm.h
View file
Open in desktop
Expand Up @@ -5,80 +5,176 @@
//
// Helper Class for applying regression-based energy corrections with optimized BDT implementation
//
// Authors: J.Bendavid
// Original Author: J.Bendavid
//
// Refactored, modernised and extended to HGCAL by S. Harper (RAL/CERN)
// with input from S. Bhattacharya (DESY)
//--------------------------------------------------------------------------------------------------

#ifndef RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorSemiParm_h
#define RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorSemiParm_h

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

#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "CondFormats/GBRForest/interface/GBRForestD.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/ESGetToken.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "CondFormats/GBRForest/interface/GBRForestD.h"
#include "CondFormats/DataRecord/interface/GBRDWrapperRcd.h"
#include "RecoEgamma/EgammaTools/interface/EgammaBDTOutputTransformer.h"
#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h"

class SCEnergyCorrectorSemiParm {
public:
SCEnergyCorrectorSemiParm();
~SCEnergyCorrectorSemiParm();
//if you want override the default on where conditions are consumed, you need to use
//the other constructor and then call setTokens approprately
SCEnergyCorrectorSemiParm(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc);

static void fillPSetDescription(edm::ParameterSetDescription& desc);
static edm::ParameterSetDescription makePSetDescription();

template <edm::Transition tr = edm::Transition::BeginLuminosityBlock>
void setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc);

void setTokens(const edm::ParameterSet &iConfig, edm::ConsumesCollector &cc);
void setEventSetup(const edm::EventSetup& es);
void setEvent(const edm::Event& e);

std::pair<double, double> getCorrections(const reco::SuperCluster& sc) const;
void modifyObject(reco::SuperCluster& sc) const;

std::vector<float> getRegData(const reco::SuperCluster& sc) const;

private:
class RegParam {
public:
RegParam(std::string meanKey = "",
float meanLow = 0,
Comment on lines +61 to +62
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
RegParam(std::string meanKey = "",
float meanLow = 0,
RegParam(std::string meanKey,
float meanLow,

from the use cases I don't see why default values are needed

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the result of two styles clashing. I prefer an initialise at construction approach while SCEnergyCorrectorSemiPara is a construct then initialise paradigm.

RegParam contains at EgammaBDTOutputTransformer which is has no way of altering its parameters after construction. However the RegParam objects only get their parameters latter. Therefore some defaults must be sent. I could just hard code the defaults into the constructor but to me it makes very little difference and in this way I can eventually adapt the class to my prefered paradigm.

Note this is the entire reason EgammaBTDOutputTransformer lost the constness of its two member variables, I needed the move constructor active and that was preferable to writing a custom one.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I perhaps missed the line(s) in code where RegParam is used without defaults. Please clarify.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I perhaps missed the line(s) in code where RegParam is used without defaults. Please clarify.

sorry, I meant with defaults

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

line 39, the default constructor of SCEnergyCorrectorSemiParm.

So I could get rid of that. But the problem is now we have to declare our transitions for ESProducts which I achieve through templates. But I dont want to template the class, just the method which calls it.

Now there are other solutions but are a few default parameters really so bad? I'm not seeing the problem.

float meanHigh = 0,
std::string sigmaKey = "",
float sigmaLow = 0,
float sigmaHigh = 0)
: meanKey_(std::move(meanKey)),
sigmaKey_(std::move(sigmaKey)),
meanOutTrans_(meanLow, meanHigh),
sigmaOutTrans_(sigmaLow, sigmaHigh) {}
RegParam(edm::ConsumesCollector cc,
std::string meanKey = "",
float meanLow = 0,
float meanHigh = 0,
std::string sigmaKey = "",
float sigmaLow = 0,
float sigmaHigh = 0)
: RegParam(meanKey, meanLow, meanHigh, sigmaKey, sigmaLow, sigmaHigh) {
setTokens(cc);
}
template <edm::Transition esTransition = edm::Transition::BeginLuminosityBlock>
void setTokens(edm::ConsumesCollector cc);
void setForests(const edm::EventSetup& setup);

std::pair<double, double> getCorrections(const reco::SuperCluster &sc) const;
void modifyObject(reco::SuperCluster &sc);
double mean(const std::vector<float>& data) const;
double sigma(const std::vector<float>& data) const;

void setEventSetup(const edm::EventSetup &es);
void setEvent(const edm::Event &e);
private:
std::string meanKey_;
std::string sigmaKey_;
EgammaBDTOutputTransformer meanOutTrans_;
EgammaBDTOutputTransformer sigmaOutTrans_;
const GBRForestD* meanForest_;
const GBRForestD* sigmaForest_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> meanForestToken_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> sigmaForestToken_;
};

protected:
const GBRForestD *foresteb_;
const GBRForestD *forestee_;
const GBRForestD *forestsigmaeb_;
const GBRForestD *forestsigmaee_;
//returns barrel for ecal barrel, otherwise returns endcap
const RegParam& getRegParam(const DetId& detId) const {
return detId.det() == DetId::Ecal && detId.subdetId() == EcalBarrel ? regParamBarrel_ : regParamEndcap_;
}

edm::ESHandle<CaloTopology> calotopo_;
edm::ESHandle<CaloGeometry> calogeom_;
std::vector<float> getRegDataECALV1(const reco::SuperCluster& sc) const;
std::vector<float> getRegDataECALHLTV1(const reco::SuperCluster& sc) const;
std::vector<float> getRegDataHGCALV1(const reco::SuperCluster& sc) const;
std::vector<float> getRegDataHGCALHLTV1(const reco::SuperCluster& sc) const;

edm::ESGetToken<CaloTopology, CaloTopologyRecord> tokenCaloTopo_;
edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokenCaloGeom_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> tokenRegressionKeyEB_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> tokenUncertaintyKeyEB_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> tokenRegressionKeyEE_;
edm::ESGetToken<GBRForestD, GBRDWrapperRcd> tokenUncertaintyKeyEE_;
//barrel = always ecal barrel, endcap may be ECAL or HGCAL
RegParam regParamBarrel_;
RegParam regParamEndcap_;

const CaloTopology* caloTopo_;
const CaloGeometry* caloGeom_;
edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopoToken_;
edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;

edm::EDGetTokenT<EcalRecHitCollection> tokenEBRecHits_;
edm::EDGetTokenT<EcalRecHitCollection> tokenEERecHits_;
edm::EDGetTokenT<reco::PFRecHitCollection> tokenHgcalRecHits_;
edm::EDGetTokenT<reco::VertexCollection> tokenVertices_;

edm::Handle<EcalRecHitCollection> recHitsEB_;
edm::Handle<EcalRecHitCollection> recHitsEE_;
edm::Handle<reco::PFRecHitCollection> recHitsHgcal_;
edm::Handle<reco::VertexCollection> vertices_;
edm::Handle<EcalRecHitCollection> rechitsEB_;
edm::Handle<EcalRecHitCollection> rechitsEE_;

edm::InputTag ecalHitsEBInputTag_;
edm::InputTag ecalHitsEEInputTag_;
edm::InputTag vertexInputTag_;

std::string regressionKeyEB_;
std::string uncertaintyKeyEB_;
std::string regressionKeyEE_;
std::string uncertaintyKeyEE_;

private:
bool isHLT_;
bool isPhaseII_;
bool applySigmaIetaIphiBug_; //there was a bug in sigmaIetaIphi for the 74X application
int nHitsAboveThreshold_;
float eThreshold_;
int nHitsAboveThresholdEB_;
int nHitsAboveThresholdEE_;
int nHitsAboveThresholdHG_;
float hitsEnergyThreshold_;
float hgcalCylinderR_;
HGCalShowerShapeHelper hgcalShowerShapes_;
};

template <edm::Transition esTransition>
void SCEnergyCorrectorSemiParm::RegParam::setTokens(edm::ConsumesCollector cc) {
meanForestToken_ = cc.esConsumes<GBRForestD, GBRDWrapperRcd, esTransition>(edm::ESInputTag("", meanKey_));
sigmaForestToken_ = cc.esConsumes<GBRForestD, GBRDWrapperRcd, esTransition>(edm::ESInputTag("", sigmaKey_));
}

template <edm::Transition esTransition>
void SCEnergyCorrectorSemiParm::setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc) {
isHLT_ = iConfig.getParameter<bool>("isHLT");
isPhaseII_ = iConfig.getParameter<bool>("isPhaseII");
applySigmaIetaIphiBug_ = iConfig.getParameter<bool>("applySigmaIetaIphiBug");
tokenEBRecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRecHitsEB"));
if (not isPhaseII_) {
tokenEERecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRecHitsEE"));
} else {
tokenHgcalRecHits_ = cc.consumes<reco::PFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hgcalRecHits"));
hgcalCylinderR_ = iConfig.getParameter<double>("hgcalCylinderR");
hgcalShowerShapes_.setTokens<esTransition>(cc);
}
caloGeomToken_ = cc.esConsumes<CaloGeometry, CaloGeometryRecord, esTransition>();
caloTopoToken_ = cc.esConsumes<CaloTopology, CaloTopologyRecord, esTransition>();

regParamBarrel_ = RegParam(iConfig.getParameter<std::string>("regressionKeyEB"),
iConfig.getParameter<double>("regressionMinEB"),
iConfig.getParameter<double>("regressionMaxEB"),
iConfig.getParameter<std::string>("uncertaintyKeyEB"),
iConfig.getParameter<double>("uncertaintyMinEB"),
iConfig.getParameter<double>("uncertaintyMaxEB"));
regParamBarrel_.setTokens<esTransition>(cc);
regParamEndcap_ = RegParam(iConfig.getParameter<std::string>("regressionKeyEE"),
iConfig.getParameter<double>("regressionMinEE"),
iConfig.getParameter<double>("regressionMaxEE"),
iConfig.getParameter<std::string>("uncertaintyKeyEE"),
iConfig.getParameter<double>("uncertaintyMinEE"),
iConfig.getParameter<double>("uncertaintyMaxEE"));
regParamEndcap_.setTokens<esTransition>(cc);
hitsEnergyThreshold_ = iConfig.getParameter<double>("eRecHitThreshold");
if (not isHLT_) {
tokenVertices_ = cc.consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexCollection"));
}
}
#endif