Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:24:38

0001 //--------------------------------------------------------------------------------------------------
0002 // $Id $
0003 //
0004 // SCEnergyCorrectorSemiParm
0005 //
0006 // Helper Class for applying regression-based energy corrections with optimized BDT implementation
0007 //
0008 // Original Author: J.Bendavid
0009 //
0010 // Refactored, modernised and extended to HGCAL by S. Harper (RAL/CERN)
0011 // with input from S. Bhattacharya (DESY)
0012 //--------------------------------------------------------------------------------------------------
0013 
0014 #ifndef RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorSemiParm_h
0015 #define RecoEcal_EgammaClusterAlgos_SCEnergyCorrectorSemiParm_h
0016 
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020 #include "FWCore/Framework/interface/ConsumesCollector.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0023 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0026 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0027 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0028 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0029 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032 #include "CondFormats/GBRForest/interface/GBRForestD.h"
0033 #include "CondFormats/DataRecord/interface/GBRDWrapperRcd.h"
0034 #include "RecoEgamma/EgammaTools/interface/EgammaBDTOutputTransformer.h"
0035 #include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h"
0036 
0037 class SCEnergyCorrectorSemiParm {
0038 public:
0039   SCEnergyCorrectorSemiParm();
0040   //if you want override the default on where conditions are consumed, you need to use
0041   //the other constructor and then call setTokens approprately
0042   SCEnergyCorrectorSemiParm(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc);
0043 
0044   static void fillPSetDescription(edm::ParameterSetDescription& desc);
0045   static edm::ParameterSetDescription makePSetDescription();
0046 
0047   template <edm::Transition tr = edm::Transition::BeginLuminosityBlock>
0048   void setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc);
0049 
0050   void setEventSetup(const edm::EventSetup& es);
0051   void setEvent(const edm::Event& e);
0052 
0053   std::pair<double, double> getCorrections(const reco::SuperCluster& sc) const;
0054   void modifyObject(reco::SuperCluster& sc) const;
0055 
0056   std::vector<float> getRegData(const reco::SuperCluster& sc) const;
0057 
0058 private:
0059   class RegParam {
0060   public:
0061     RegParam(std::string meanKey = "",
0062              float meanLow = 0,
0063              float meanHigh = 0,
0064              std::string sigmaKey = "",
0065              float sigmaLow = 0,
0066              float sigmaHigh = 0)
0067         : meanKey_(std::move(meanKey)),
0068           sigmaKey_(std::move(sigmaKey)),
0069           meanOutTrans_(meanLow, meanHigh),
0070           sigmaOutTrans_(sigmaLow, sigmaHigh) {}
0071     RegParam(edm::ConsumesCollector cc,
0072              std::string meanKey = "",
0073              float meanLow = 0,
0074              float meanHigh = 0,
0075              std::string sigmaKey = "",
0076              float sigmaLow = 0,
0077              float sigmaHigh = 0)
0078         : RegParam(meanKey, meanLow, meanHigh, sigmaKey, sigmaLow, sigmaHigh) {
0079       setTokens(cc);
0080     }
0081     template <edm::Transition esTransition = edm::Transition::BeginLuminosityBlock>
0082     void setTokens(edm::ConsumesCollector cc);
0083     void setForests(const edm::EventSetup& setup);
0084 
0085     double mean(const std::vector<float>& data) const;
0086     double sigma(const std::vector<float>& data) const;
0087 
0088   private:
0089     std::string meanKey_;
0090     std::string sigmaKey_;
0091     EgammaBDTOutputTransformer meanOutTrans_;
0092     EgammaBDTOutputTransformer sigmaOutTrans_;
0093     const GBRForestD* meanForest_;
0094     const GBRForestD* sigmaForest_;
0095     edm::ESGetToken<GBRForestD, GBRDWrapperRcd> meanForestToken_;
0096     edm::ESGetToken<GBRForestD, GBRDWrapperRcd> sigmaForestToken_;
0097   };
0098 
0099   //returns barrel for ecal barrel, otherwise returns endcap
0100   const RegParam& getRegParam(const DetId& detId) const {
0101     return detId.det() == DetId::Ecal && detId.subdetId() == EcalBarrel ? regParamBarrel_ : regParamEndcap_;
0102   }
0103 
0104   std::vector<float> getRegDataECALV1(const reco::SuperCluster& sc) const;
0105   std::vector<float> getRegDataECALHLTV1(const reco::SuperCluster& sc) const;
0106   std::vector<float> getRegDataHGCALV1(const reco::SuperCluster& sc) const;
0107   std::vector<float> getRegDataHGCALHLTV1(const reco::SuperCluster& sc) const;
0108 
0109   //barrel = always ecal barrel, endcap may be ECAL or HGCAL
0110   RegParam regParamBarrel_;
0111   RegParam regParamEndcap_;
0112 
0113   const CaloTopology* caloTopo_;
0114   const CaloGeometry* caloGeom_;
0115   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopoToken_;
0116   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0117 
0118   edm::EDGetTokenT<EcalRecHitCollection> tokenEBRecHits_;
0119   edm::EDGetTokenT<EcalRecHitCollection> tokenEERecHits_;
0120   edm::EDGetTokenT<reco::PFRecHitCollection> tokenHgcalRecHits_;
0121   edm::EDGetTokenT<reco::VertexCollection> tokenVertices_;
0122 
0123   edm::Handle<EcalRecHitCollection> recHitsEB_;
0124   edm::Handle<EcalRecHitCollection> recHitsEE_;
0125   edm::Handle<reco::PFRecHitCollection> recHitsHgcal_;
0126   edm::Handle<reco::VertexCollection> vertices_;
0127 
0128   bool isHLT_;
0129   bool isPhaseII_;
0130   bool regTrainedWithPS_;
0131   bool applySigmaIetaIphiBug_;  //there was a bug in sigmaIetaIphi for the 74X application
0132   int nHitsAboveThresholdEB_;
0133   int nHitsAboveThresholdEE_;
0134   int nHitsAboveThresholdHG_;
0135   float hitsEnergyThreshold_;
0136   float hgcalCylinderR_;
0137   HGCalShowerShapeHelper hgcalShowerShapes_;
0138 };
0139 
0140 template <edm::Transition esTransition>
0141 void SCEnergyCorrectorSemiParm::RegParam::setTokens(edm::ConsumesCollector cc) {
0142   meanForestToken_ = cc.esConsumes<GBRForestD, GBRDWrapperRcd, esTransition>(edm::ESInputTag("", meanKey_));
0143   sigmaForestToken_ = cc.esConsumes<GBRForestD, GBRDWrapperRcd, esTransition>(edm::ESInputTag("", sigmaKey_));
0144 }
0145 
0146 template <edm::Transition esTransition>
0147 void SCEnergyCorrectorSemiParm::setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector cc) {
0148   isHLT_ = iConfig.getParameter<bool>("isHLT");
0149   isPhaseII_ = iConfig.getParameter<bool>("isPhaseII");
0150   regTrainedWithPS_ = iConfig.getParameter<bool>("regTrainedWithPS");
0151   applySigmaIetaIphiBug_ = iConfig.getParameter<bool>("applySigmaIetaIphiBug");
0152   tokenEBRecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRecHitsEB"));
0153   if (not isPhaseII_) {
0154     tokenEERecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalRecHitsEE"));
0155   } else {
0156     tokenHgcalRecHits_ = cc.consumes<reco::PFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hgcalRecHits"));
0157     hgcalCylinderR_ = iConfig.getParameter<double>("hgcalCylinderR");
0158     hgcalShowerShapes_.setTokens<esTransition>(cc);
0159   }
0160   caloGeomToken_ = cc.esConsumes<CaloGeometry, CaloGeometryRecord, esTransition>();
0161   caloTopoToken_ = cc.esConsumes<CaloTopology, CaloTopologyRecord, esTransition>();
0162 
0163   regParamBarrel_ = RegParam(iConfig.getParameter<std::string>("regressionKeyEB"),
0164                              iConfig.getParameter<double>("regressionMinEB"),
0165                              iConfig.getParameter<double>("regressionMaxEB"),
0166                              iConfig.getParameter<std::string>("uncertaintyKeyEB"),
0167                              iConfig.getParameter<double>("uncertaintyMinEB"),
0168                              iConfig.getParameter<double>("uncertaintyMaxEB"));
0169   regParamBarrel_.setTokens<esTransition>(cc);
0170   regParamEndcap_ = RegParam(iConfig.getParameter<std::string>("regressionKeyEE"),
0171                              iConfig.getParameter<double>("regressionMinEE"),
0172                              iConfig.getParameter<double>("regressionMaxEE"),
0173                              iConfig.getParameter<std::string>("uncertaintyKeyEE"),
0174                              iConfig.getParameter<double>("uncertaintyMinEE"),
0175                              iConfig.getParameter<double>("uncertaintyMaxEE"));
0176   regParamEndcap_.setTokens<esTransition>(cc);
0177   hitsEnergyThreshold_ = iConfig.getParameter<double>("eRecHitThreshold");
0178   if (not isHLT_) {
0179     tokenVertices_ = cc.consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexCollection"));
0180   }
0181 }
0182 #endif