Go to the documentation of this file.00001 #include "RecoEcal/EgammaClusterProducers/interface/EgammaSCCorrectionMaker.h"
00002 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00003 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00004 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00006
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 #include "FWCore/Utilities/interface/Exception.h"
00009 #include "DataFormats/Common/interface/Handle.h"
00010
00011 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00013 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00014 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00015 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00016 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
00017 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
00018
00019 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
00020 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
00021
00022 #include <string>
00023
00024 EgammaSCCorrectionMaker::EgammaSCCorrectionMaker(const edm::ParameterSet& ps)
00025 {
00026
00027
00028 std::string debugString = ps.getParameter<std::string>("VerbosityLevel");
00029 if (debugString == "DEBUG") verbosity_ = EgammaSCEnergyCorrectionAlgo::pDEBUG;
00030 else if (debugString == "INFO") verbosity_ = EgammaSCEnergyCorrectionAlgo::pINFO;
00031 else verbosity_ = EgammaSCEnergyCorrectionAlgo::pERROR;
00032
00033
00034 rHInputProducer_ = ps.getParameter<edm::InputTag>("recHitProducer");
00035 sCInputProducer_ = ps.getParameter<edm::InputTag>("rawSuperClusterProducer");
00036 std::string sCAlgo_str = ps.getParameter<std::string>("superClusterAlgo");
00037
00038
00039
00040 edm::ParameterSet fCorrPset;
00041 if (sCAlgo_str=="Hybrid") {
00042 sCAlgo_= reco::CaloCluster::hybrid;
00043 fCorrPset = ps.getParameter<edm::ParameterSet>("hyb_fCorrPset");
00044 } else if (sCAlgo_str=="Island") {
00045 sCAlgo_= reco::CaloCluster::island;
00046 fCorrPset = ps.getParameter<edm::ParameterSet>("isl_fCorrPset");
00047 } else if (sCAlgo_str=="DynamicHybrid") {
00048 sCAlgo_ = reco::CaloCluster::dynamicHybrid;
00049 fCorrPset = ps.getParameter<edm::ParameterSet>("dyn_fCorrPset");
00050 } else if (sCAlgo_str=="Multi5x5") {
00051 sCAlgo_ = reco::CaloCluster::multi5x5;
00052 fCorrPset = ps.getParameter<edm::ParameterSet>("fix_fCorrPset");
00053 } else {
00054 edm::LogError("EgammaSCCorrectionMakerError")
00055 << "Error! SuperClusterAlgo in config file must be Hybrid or Island: "
00056 << sCAlgo_str << " Using Hybrid by default";
00057 sCAlgo_=reco::CaloCluster::hybrid;
00058 }
00059
00060
00061 applyEnergyCorrection_ = ps.getParameter<bool>("applyEnergyCorrection");
00062 applyCrackCorrection_ = ps.getParameter<bool>("applyCrackCorrection");
00063
00064 sigmaElectronicNoise_ = ps.getParameter<double>("sigmaElectronicNoise");
00065
00066 etThresh_ = ps.getParameter<double>("etThresh");
00067
00068
00069 outputCollection_ = ps.getParameter<std::string>("corectedSuperClusterCollection");
00070 produces<reco::SuperClusterCollection>(outputCollection_);
00071
00072
00073 energyCorrector_ = new EgammaSCEnergyCorrectionAlgo(sigmaElectronicNoise_, sCAlgo_, fCorrPset, verbosity_);
00074
00075
00076 if (applyEnergyCorrection_ )
00077 energyCorrectionFunction_ = EcalClusterFunctionFactory::get()->create("EcalClusterEnergyCorrection", ps);
00078 else
00079 energyCorrectionFunction_=0;
00080 if (applyCrackCorrection_ )
00081 crackCorrectionFunction_ = EcalClusterFunctionFactory::get()->create("EcalClusterCrackCorrection", ps);
00082 else
00083 crackCorrectionFunction_=0;
00084
00085
00086 }
00087
00088 EgammaSCCorrectionMaker::~EgammaSCCorrectionMaker()
00089 {
00090 if (energyCorrectionFunction_) delete energyCorrectionFunction_;
00091 if (crackCorrectionFunction_) delete crackCorrectionFunction_;
00092 delete energyCorrector_;
00093 }
00094
00095 void
00096 EgammaSCCorrectionMaker::produce(edm::Event& evt, const edm::EventSetup& es)
00097 {
00098 using namespace edm;
00099
00100
00101 if(applyEnergyCorrection_)
00102 energyCorrectionFunction_->init(es);
00103
00104
00105 if(applyCrackCorrection_)
00106 crackCorrectionFunction_->init(es);
00107
00108
00109
00110 edm::ESHandle<CaloGeometry> geoHandle;
00111 es.get<CaloGeometryRecord>().get(geoHandle);
00112 const CaloGeometry& geometry = *geoHandle;
00113 const CaloSubdetectorGeometry *geometry_p;
00114
00115 std::string rHInputCollection = rHInputProducer_.instance();
00116 if(rHInputCollection == "EcalRecHitsEB") {
00117 geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00118 } else if(rHInputCollection == "EcalRecHitsEE") {
00119 geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00120 } else if(rHInputCollection == "EcalRecHitsPS") {
00121 geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
00122 } else {
00123 std::string str = "\n\nSCCorrectionMaker encountered invalied ecalhitcollection type: " + rHInputCollection + ".\n\n";
00124 throw(std::runtime_error( str.c_str() ));
00125 }
00126
00127
00128 Handle<reco::SuperClusterCollection> pRawSuperClusters;
00129 try {
00130 evt.getByLabel(sCInputProducer_, pRawSuperClusters);
00131 } catch ( cms::Exception& ex ) {
00132 edm::LogError("EgammaSCCorrectionMakerError")
00133 << "Error! can't get the rawSuperClusters "
00134 << sCInputProducer_.label() ;
00135 }
00136
00137
00138 Handle<EcalRecHitCollection> pRecHits;
00139 try {
00140 evt.getByLabel(rHInputProducer_, pRecHits);
00141 } catch ( cms::Exception& ex ) {
00142 edm::LogError("EgammaSCCorrectionMakerError")
00143 << "Error! can't get the RecHits "
00144 << rHInputProducer_.label();
00145 }
00146
00147
00148 const EcalRecHitCollection *hitCollection = pRecHits.product();
00149 const reco::SuperClusterCollection *rawClusters = pRawSuperClusters.product();
00150
00151
00152 std::auto_ptr<reco::SuperClusterCollection> corrClusters(new reco::SuperClusterCollection);
00153
00154
00155 reco::SuperClusterCollection::const_iterator aClus;
00156 for(aClus = rawClusters->begin(); aClus != rawClusters->end(); aClus++)
00157 {
00158 reco::SuperCluster enecorrClus,crackcorrClus;
00159
00160 if(applyEnergyCorrection_)
00161 enecorrClus = energyCorrector_->applyCorrection(*aClus, *hitCollection, sCAlgo_, geometry_p, energyCorrectionFunction_);
00162 else
00163 enecorrClus=*aClus;
00164
00165 if(applyCrackCorrection_)
00166 crackcorrClus=energyCorrector_->applyCrackCorrection(enecorrClus,crackCorrectionFunction_);
00167 else
00168 crackcorrClus=enecorrClus;
00169
00170 if(crackcorrClus.energy()*sin(crackcorrClus.position().theta())>etThresh_) {
00171
00172 corrClusters->push_back(crackcorrClus);
00173 }
00174 }
00175
00176 evt.put(corrClusters, outputCollection_);
00177
00178 }
00179