ElectronTools/plugins/ElectronRegressionEnergyProducer.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/Framework/interface/stream/EDFilter.h"
0015
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0018 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0019 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0020 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0021 #include "DataFormats/PatCandidates/interface/Electron.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyRegressionEvaluate.h"
0025
0026 //
0027 // class declaration
0028 //
0029
0030 class ElectronRegressionEnergyProducer : public edm::stream::EDFilter<> {
0031 public:
0032   explicit ElectronRegressionEnergyProducer(const edm::ParameterSet&);
0033   ~ElectronRegressionEnergyProducer() override;
0034
0035 private:
0036   bool filter(edm::Event&, const edm::EventSetup&) override;
0037
0038   // ----------member data ---------------------------
0039   bool printDebug_;
0040   edm::EDGetTokenT<reco::GsfElectronCollection> electronToken_;
0041
0042   std::string regressionInputFile_;
0043   uint32_t energyRegressionType_;
0044
0045   std::string nameEnergyReg_;
0046   std::string nameEnergyErrorReg_;
0047
0048   bool geomInitialized_;
0049
0050   const CaloTopology* ecalTopology_;
0051   const CaloGeometry* caloGeometry_;
0052
0053   edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEBToken_;
0054   edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEEToken_;
0055
0056   ElectronEnergyRegressionEvaluate* regressionEvaluator;
0057
0058   edm::EDGetTokenT<reco::VertexCollection> hVertexToken_;
0059   edm::EDGetTokenT<double> hRhoKt6PFJetsToken_;
0060
0061   edm::ESGetToken<CaloTopology, CaloTopologyRecord> ecalTopoToken_;
0062   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0063 };
0064
0065 ElectronRegressionEnergyProducer::ElectronRegressionEnergyProducer(const edm::ParameterSet& iConfig) {
0066   printDebug_ = iConfig.getUntrackedParameter<bool>("printDebug", false);
0067   electronToken_ = consumes<reco::GsfElectronCollection>(iConfig.getParameter<edm::InputTag>("electronTag"));
0068
0069   regressionInputFile_ = iConfig.getParameter<std::string>("regressionInputFile");
0070   energyRegressionType_ = iConfig.getParameter<uint32_t>("energyRegressionType");
0071
0072   nameEnergyReg_ = iConfig.getParameter<std::string>("nameEnergyReg");
0073   nameEnergyErrorReg_ = iConfig.getParameter<std::string>("nameEnergyErrorReg");
0074
0075   recHitCollectionEBToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitCollectionEB"));
0076   recHitCollectionEEToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitCollectionEE"));
0077
0078   hVertexToken_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0079   hRhoKt6PFJetsToken_ = consumes<double>(edm::InputTag("kt6PFJets", "rho"));
0080
0081   ecalTopoToken_ = esConsumes();
0082   caloGeomToken_ = esConsumes();
0083
0084   produces<edm::ValueMap<double> >(nameEnergyReg_);
0085   produces<edm::ValueMap<double> >(nameEnergyErrorReg_);
0086
0087   regressionEvaluator = new ElectronEnergyRegressionEvaluate();
0088
0089   //set regression type
0090   ElectronEnergyRegressionEvaluate::ElectronEnergyRegressionType type = ElectronEnergyRegressionEvaluate::kNoTrkVar;
0091   if (energyRegressionType_ == 1)
0092     type = ElectronEnergyRegressionEvaluate::kNoTrkVar;
0093   else if (energyRegressionType_ == 2)
0094     type = ElectronEnergyRegressionEvaluate::kWithSubCluVar;
0095   else if (energyRegressionType_ == 3)
0096     type = ElectronEnergyRegressionEvaluate::kWithTrkVarV1;
0097   else if (energyRegressionType_ == 4)
0098     type = ElectronEnergyRegressionEvaluate::kWithTrkVarV2;
0099
0100   //load weights and initialize
0101   regressionEvaluator->initialize(regressionInputFile_, type);
0102
0103   geomInitialized_ = false;
0104 }
0105
0106 ElectronRegressionEnergyProducer::~ElectronRegressionEnergyProducer() { delete regressionEvaluator; }
0107
0108 // ------------ method called on each new Event  ------------
0109 bool ElectronRegressionEnergyProducer::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110   assert(regressionEvaluator->isInitialized());
0111
0112   if (!geomInitialized_) {
0113     ecalTopology_ = &(iSetup.getData(ecalTopoToken_));
0114     caloGeometry_ = &(iSetup.getData(caloGeomToken_));
0115     geomInitialized_ = true;
0116   }
0117
0118   std::unique_ptr<edm::ValueMap<double> > regrEnergyMap(new edm::ValueMap<double>());
0119   edm::ValueMap<double>::Filler energyFiller(*regrEnergyMap);
0120
0121   std::unique_ptr<edm::ValueMap<double> > regrEnergyErrorMap(new edm::ValueMap<double>());
0122   edm::ValueMap<double>::Filler energyErrorFiller(*regrEnergyErrorMap);
0123
0124   edm::Handle<reco::GsfElectronCollection> egCollection;
0125   iEvent.getByToken(electronToken_, egCollection);
0126   const reco::GsfElectronCollection egCandidates = (*egCollection.product());
0127
0128   std::vector<double> energyValues;
0129   std::vector<double> energyErrorValues;
0130   energyValues.reserve(egCollection->size());
0131   energyErrorValues.reserve(egCollection->size());
0132   //
0133   //**************************************************************************
0134   // Rechits
0135   //**************************************************************************
0136   edm::Handle<EcalRecHitCollection> pEBRecHits;
0137   edm::Handle<EcalRecHitCollection> pEERecHits;
0138   iEvent.getByToken(recHitCollectionEBToken_, pEBRecHits);
0139   iEvent.getByToken(recHitCollectionEEToken_, pEERecHits);
0140
0141   //**************************************************************************
0142   //Get Number of Vertices
0143   //**************************************************************************
0144   edm::Handle<reco::VertexCollection> hVertexProduct;
0145   iEvent.getByToken(hVertexToken_, hVertexProduct);
0146   const reco::VertexCollection inVertices = *(hVertexProduct.product());
0147
0148   // loop through all vertices
0149   Int_t nvertices = 0;
0150   for (reco::VertexCollection::const_iterator inV = inVertices.begin(); inV != inVertices.end(); ++inV) {
0151     // pass these vertex cuts
0152     if (inV->ndof() >= 4 && inV->position().Rho() <= 2.0 && fabs(inV->z()) <= 24.0) {
0153       nvertices++;
0154     }
0155   }
0156
0157   //**************************************************************************
0158   //Get Rho
0159   //**************************************************************************
0160   double rho = 0;
0161   edm::Handle<double> hRhoKt6PFJets;
0162   iEvent.getByToken(hRhoKt6PFJetsToken_, hRhoKt6PFJets);
0163   rho = (*hRhoKt6PFJets);
0164
0165   for (reco::GsfElectronCollection::const_iterator egIter = egCandidates.begin(); egIter != egCandidates.end();
0166        ++egIter) {
0167     const EcalRecHitCollection* recHits = nullptr;
0168     if (egIter->isEB())
0169       recHits = pEBRecHits.product();
0170     else
0171       recHits = pEERecHits.product();
0172
0173     SuperClusterHelper mySCHelper(&(*egIter), recHits, ecalTopology_, caloGeometry_);
0174
0175     double energy = regressionEvaluator->calculateRegressionEnergy(&(*egIter), mySCHelper, rho, nvertices, printDebug_);
0176
0177     double error =
0178         regressionEvaluator->calculateRegressionEnergyUncertainty(&(*egIter), mySCHelper, rho, nvertices, printDebug_);
0179
0180     energyValues.push_back(energy);
0181     energyErrorValues.push_back(error);
0182   }
0183
0184   energyFiller.insert(egCollection, energyValues.begin(), energyValues.end());
0185   energyFiller.fill();
0186
0187   energyErrorFiller.insert(egCollection, energyErrorValues.begin(), energyErrorValues.end());
0188   energyErrorFiller.fill();
0189
0190   iEvent.put(std::move(regrEnergyMap), nameEnergyReg_);
0191   iEvent.put(std::move(regrEnergyErrorMap), nameEnergyErrorReg_);
0192
0193   return true;
0194 }
0195
0196 //define this as a plug-in
0197 DEFINE_FWK_MODULE(ElectronRegressionEnergyProducer);