Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-05 07:05:21

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);