EgammaClusterProducers/src/EgammaSCCorrectionMaker.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaSCCorrectionMaker
0004 // Class:      EgammaSCCorrectionMaker
0005 //
0006 /**\class EgammaSCCorrectionMaker EgammaSCCorrectionMaker.cc EgammaSCCorrectionMaker/EgammaSCCorrectionMaker/src/EgammaSCCorrectionMaker.cc
0007
0008  Description: Producer of corrected SuperClusters
0009
0010 */
0011 //
0012 // Original Author:  Dave Evans
0013 //         Created:  Thu Apr 13 15:50:17 CEST 2006
0014 //
0015 //
0016
0017 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0020 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0021 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0022 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0023 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Utilities/interface/Exception.h"
0030 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0032 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0033 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0034 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0035 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0036 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0037 #include "RecoEcal/EgammaClusterAlgos/interface/EgammaSCEnergyCorrectionAlgo.h"
0038 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
0039 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0040
0041 #include "EcalBasicClusterLocalContCorrection.h"
0042
0043 #include <memory>
0044 #include <string>
0045
0046 class EgammaSCCorrectionMaker : public edm::stream::EDProducer<> {
0047 public:
0048   explicit EgammaSCCorrectionMaker(const edm::ParameterSet&);
0049   void produce(edm::Event&, const edm::EventSetup&) override;
0050
0051 private:
0052   std::unique_ptr<EcalClusterFunctionBaseClass> energyCorrectionFunction_;
0053   std::unique_ptr<EcalClusterFunctionBaseClass> crackCorrectionFunction_;
0054   std::unique_ptr<EcalBasicClusterLocalContCorrection> localContCorrectionFunction_;
0055
0056   // pointer to the correction algo object
0057   std::unique_ptr<EgammaSCEnergyCorrectionAlgo> energyCorrector_;
0058
0059   // vars for the correction algo
0060   bool applyEnergyCorrection_;
0061   bool applyCrackCorrection_;
0062   bool applyLocalContCorrection_;
0063
0064   std::string energyCorrectorName_;
0065   std::string crackCorrectorName_;
0066
0067   int modeEB_;
0068   int modeEE_;
0069
0070   //     bool oldEnergyScaleCorrection_;
0071   double sigmaElectronicNoise_;
0072   double etThresh_;
0073
0074   // vars to get products
0075   edm::EDGetTokenT<EcalRecHitCollection> rHInputProducer_;
0076   edm::EDGetTokenT<reco::SuperClusterCollection> sCInputProducer_;
0077   edm::InputTag rHTag_;
0078   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0079
0080   reco::CaloCluster::AlgoId sCAlgo_;
0081   std::string outputCollection_;
0082 };
0083
0084 #include "FWCore/Framework/interface/MakerMacros.h"
0085 DEFINE_FWK_MODULE(EgammaSCCorrectionMaker);
0086
0087 EgammaSCCorrectionMaker::EgammaSCCorrectionMaker(const edm::ParameterSet& ps) {
0088   // the input producers
0089   rHTag_ = ps.getParameter<edm::InputTag>("recHitProducer");
0090   rHInputProducer_ = consumes<EcalRecHitCollection>(rHTag_);
0091   sCInputProducer_ = consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("rawSuperClusterProducer"));
0092   caloGeomToken_ = esConsumes();
0093   std::string sCAlgo_str = ps.getParameter<std::string>("superClusterAlgo");
0094
0095   // determine which BasicCluster algo we are correcting for
0096   //And obtain forrection parameters form cfg file
0097   edm::ParameterSet fCorrPset;
0098   if (sCAlgo_str == "Hybrid") {
0099     sCAlgo_ = reco::CaloCluster::hybrid;
0100     fCorrPset = ps.getParameter<edm::ParameterSet>("hyb_fCorrPset");
0101   } else if (sCAlgo_str == "Island") {
0102     sCAlgo_ = reco::CaloCluster::island;
0103     fCorrPset = ps.getParameter<edm::ParameterSet>("isl_fCorrPset");
0104   } else if (sCAlgo_str == "DynamicHybrid") {
0105     sCAlgo_ = reco::CaloCluster::dynamicHybrid;
0106     fCorrPset = ps.getParameter<edm::ParameterSet>("dyn_fCorrPset");
0107   } else if (sCAlgo_str == "Multi5x5") {
0108     sCAlgo_ = reco::CaloCluster::multi5x5;
0109     fCorrPset = ps.getParameter<edm::ParameterSet>("fix_fCorrPset");
0110   } else {
0111     edm::LogError("EgammaSCCorrectionMakerError")
0112         << "Error! SuperClusterAlgo in config file must be Hybrid or Island: " << sCAlgo_str
0113         << "  Using Hybrid by default";
0114     sCAlgo_ = reco::CaloCluster::hybrid;
0115   }
0116
0117   // set correction algo parameters
0118   applyEnergyCorrection_ = ps.getParameter<bool>("applyEnergyCorrection");
0119   applyCrackCorrection_ = ps.getParameter<bool>("applyCrackCorrection");
0120   applyLocalContCorrection_ =
0121       ps.existsAs<bool>("applyLocalContCorrection") ? ps.getParameter<bool>("applyLocalContCorrection") : false;
0122
0123   energyCorrectorName_ = ps.getParameter<std::string>("energyCorrectorName");
0124   crackCorrectorName_ = ps.existsAs<std::string>("crackCorrectorName")
0125                             ? ps.getParameter<std::string>("crackCorrectorName")
0126                             : std::string("EcalClusterCrackCorrection");
0127
0128   modeEB_ = ps.getParameter<int>("modeEB");
0129   modeEE_ = ps.getParameter<int>("modeEE");
0130
0131   sigmaElectronicNoise_ = ps.getParameter<double>("sigmaElectronicNoise");
0132
0133   etThresh_ = ps.getParameter<double>("etThresh");
0134
0135   // set the producer parameters
0136   outputCollection_ = ps.getParameter<std::string>("corectedSuperClusterCollection");
0137   produces<reco::SuperClusterCollection>(outputCollection_);
0138
0139   // instanciate the correction algo object
0140   energyCorrector_ = std::make_unique<EgammaSCEnergyCorrectionAlgo>(sigmaElectronicNoise_);
0141
0142   // energy correction class
0143   if (applyEnergyCorrection_)
0144     energyCorrectionFunction_ =
0145         EcalClusterFunctionFactory::get()->create(energyCorrectorName_, ps, consumesCollector());
0146   //energyCorrectionFunction_ = EcalClusterFunctionFactory::get()->create("EcalClusterEnergyCorrection", ps);
0147
0148   if (applyCrackCorrection_)
0149     crackCorrectionFunction_ = EcalClusterFunctionFactory::get()->create(crackCorrectorName_, ps, consumesCollector());
0150
0151   if (applyLocalContCorrection_)
0152     localContCorrectionFunction_ = std::make_unique<EcalBasicClusterLocalContCorrection>(consumesCollector());
0153 }
0154
0155 void EgammaSCCorrectionMaker::produce(edm::Event& evt, const edm::EventSetup& es) {
0156   using namespace edm;
0157
0158   // initialize energy correction class
0159   if (applyEnergyCorrection_)
0160     energyCorrectionFunction_->init(es);
0161
0162   // initialize energy correction class
0163   if (applyCrackCorrection_)
0164     crackCorrectionFunction_->init(es);
0165
0166   // initialize containemnt correction class
0167   if (applyLocalContCorrection_)
0168     localContCorrectionFunction_->init(es);
0169
0170   // get the collection geometry:
0171   const CaloGeometry& geometry = es.getData(caloGeomToken_);
0172   const CaloSubdetectorGeometry* geometry_p;
0173
0174   std::string rHInputCollection = rHTag_.instance();
0175   if (rHInputCollection == "EcalRecHitsEB") {
0176     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0177   } else if (rHInputCollection == "EcalRecHitsEE") {
0178     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0179   } else if (rHInputCollection == "EcalRecHitsPS") {
0180     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0181   } else {
0182     std::string str =
0183         "\n\nSCCorrectionMaker encountered invalied ecalhitcollection type: " + rHInputCollection + ".\n\n";
0184     throw(std::runtime_error(str.c_str()));
0185   }
0186
0187   // Get raw SuperClusters from the event
0188   Handle<reco::SuperClusterCollection> pRawSuperClusters;
0189   evt.getByToken(sCInputProducer_, pRawSuperClusters);
0190
0191   // Get the RecHits from the event
0192   Handle<EcalRecHitCollection> pRecHits;
0193   evt.getByToken(rHInputProducer_, pRecHits);
0194
0195   // Create a pointer to the RecHits and raw SuperClusters
0196   const EcalRecHitCollection* hitCollection = pRecHits.product();
0197   const reco::SuperClusterCollection* rawClusters = pRawSuperClusters.product();
0198
0199   // Define a collection of corrected SuperClusters to put back into the event
0200   auto corrClusters = std::make_unique<reco::SuperClusterCollection>();
0201
0202   //  Loop over raw clusters and make corrected ones
0203   reco::SuperClusterCollection::const_iterator aClus;
0204   for (aClus = rawClusters->begin(); aClus != rawClusters->end(); aClus++) {
0205     reco::SuperCluster enecorrClus, crackcorrClus, localContCorrClus;
0206
0207     if (applyEnergyCorrection_)
0208       enecorrClus = energyCorrector_->applyCorrection(*aClus,
0209                                                       *hitCollection,
0210                                                       sCAlgo_,
0211                                                       geometry_p,
0212                                                       energyCorrectionFunction_.get(),
0213                                                       energyCorrectorName_,
0214                                                       modeEB_,
0215                                                       modeEE_);
0216     else
0217       enecorrClus = *aClus;
0218
0219     if (applyCrackCorrection_)
0220       crackcorrClus = EgammaSCEnergyCorrectionAlgo::applyCrackCorrection(enecorrClus, crackCorrectionFunction_.get());
0221     else
0222       crackcorrClus = enecorrClus;
0223
0224     if (applyLocalContCorrection_)
0225       localContCorrClus =
0226           EgammaSCEnergyCorrectionAlgo::applyLocalContCorrection(crackcorrClus, *localContCorrectionFunction_);
0227     else
0228       localContCorrClus = crackcorrClus;
0229
0230     if (localContCorrClus.energy() * sin(localContCorrClus.position().theta()) > etThresh_) {
0231       corrClusters->push_back(localContCorrClus);
0232     }
0233   }
0234
0235   // Put collection of corrected SuperClusters into the event
0236   evt.put(std::move(corrClusters), outputCollection_);
0237 }