File indexing completed on 2024-04-06 12:24:58
0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/PhotonEnergyCorrector.h"
0002 #include "RecoEgamma/EgammaTools/interface/egEnergyCorrectorFactoryFromRootFile.h"
0003 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0004 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0007 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0008 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0009 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011
0012 #include "RecoEgamma/EgammaPhotonAlgos/interface/EnergyUncertaintyPhotonSpecific.h"
0013
0014 PhotonEnergyCorrector::PhotonEnergyCorrector(const edm::ParameterSet& config, edm::ConsumesCollector&& iC)
0015 : ecalClusterToolsESGetTokens_{iC}, caloGeomToken_{iC.esConsumes()} {
0016 minR9Barrel_ = config.getParameter<double>("minR9Barrel");
0017 minR9Endcap_ = config.getParameter<double>("minR9Endcap");
0018
0019
0020 barrelEcalHits_ = config.getParameter<edm::InputTag>("barrelEcalHits");
0021 endcapEcalHits_ = config.getParameter<edm::InputTag>("endcapEcalHits");
0022 barrelEcalHitsToken_ = iC.consumes<EcalRecHitCollection>(config.getParameter<edm::InputTag>("barrelEcalHits"));
0023 endcapEcalHitsToken_ = iC.consumes<EcalRecHitCollection>(config.getParameter<edm::InputTag>("endcapEcalHits"));
0024
0025
0026
0027
0028 std::string superClusterFunctionName = config.getParameter<std::string>("superClusterEnergyCorrFunction");
0029 scEnergyFunction_ = EcalClusterFunctionFactory::get()->create(superClusterFunctionName, config, iC);
0030
0031
0032 std::string superClusterCrackFunctionName = config.getParameter<std::string>("superClusterCrackEnergyCorrFunction");
0033 scCrackEnergyFunction_ = EcalClusterFunctionFactory::get()->create(superClusterCrackFunctionName, config, iC);
0034
0035
0036 std::string superClusterErrorFunctionName = config.getParameter<std::string>("superClusterEnergyErrorFunction");
0037 scEnergyErrorFunction_ = EcalClusterFunctionFactory::get()->create(superClusterErrorFunctionName, config, iC);
0038
0039
0040 std::string photonEnergyFunctionName = config.getParameter<std::string>("photonEcalEnergyCorrFunction");
0041 photonEcalEnergyCorrFunction_ = EcalClusterFunctionFactory::get()->create(photonEnergyFunctionName, config, iC);
0042
0043 photonUncertaintyCalculator_ = std::make_unique<EnergyUncertaintyPhotonSpecific>(config);
0044
0045 if (config.existsAs<edm::ParameterSet>("regressionConfig")) {
0046 const edm::ParameterSet& regr_conf = config.getParameterSet("regressionConfig");
0047 const std::string& mname = regr_conf.getParameter<std::string>("modifierName");
0048 gedRegression_ = ModifyObjectValueFactory::get()->create(mname, regr_conf, iC);
0049 }
0050
0051
0052 weightsfromDB_ = config.getParameter<bool>("regressionWeightsFromDB");
0053 w_file_ = config.getParameter<std::string>("energyRegressionWeightsFileLocation");
0054 if (weightsfromDB_) {
0055 w_db_ = config.getParameter<std::string>("energyRegressionWeightsDBLocation");
0056 regressionCorrectorFactory_ = std::make_unique<EGEnergyCorrectorFactoryFromEventSetup>(iC, w_db_);
0057 } else if (w_file_ != "none") {
0058 regressionCorrector_ = std::make_unique<EGEnergyCorrector>(egEnergyCorrectorFactoryFromRootFile(w_file_.c_str()));
0059 } else {
0060 regressionCorrector_ = std::make_unique<EGEnergyCorrector>();
0061 }
0062 }
0063
0064 void PhotonEnergyCorrector::init(const edm::EventSetup& theEventSetup) {
0065 theCaloGeom_ = theEventSetup.getHandle(caloGeomToken_);
0066
0067 scEnergyFunction_->init(theEventSetup);
0068 scCrackEnergyFunction_->init(theEventSetup);
0069 scEnergyErrorFunction_->init(theEventSetup);
0070 photonEcalEnergyCorrFunction_->init(theEventSetup);
0071
0072 if (not regressionCorrector_) {
0073 regressionCorrector_ = std::make_unique<EGEnergyCorrector>(regressionCorrectorFactory_->build(theEventSetup));
0074 }
0075 photonUncertaintyCalculator_->init(theEventSetup);
0076 }
0077
0078 void PhotonEnergyCorrector::calculate(edm::Event& evt,
0079 reco::Photon& thePhoton,
0080 int subdet,
0081 const reco::VertexCollection& vtxcol,
0082 const edm::EventSetup& iSetup) {
0083 double phoEcalEnergy = -9999.;
0084 double phoEcalEnergyError = -9999.;
0085 double phoRegr1Energy = -9999.;
0086 double phoRegr1EnergyError = -9999.;
0087 theCaloGeom_->getSubdetectorGeometry(DetId::Ecal, subdet);
0088
0089 double minR9 = 0;
0090 if (subdet == EcalBarrel) {
0091 minR9 = minR9Barrel_;
0092 } else if (subdet == EcalEndcap) {
0093 minR9 = minR9Endcap_;
0094 }
0095
0096 EcalClusterLazyTools lazyTools(
0097 evt, ecalClusterToolsESGetTokens_.get(iSetup), barrelEcalHitsToken_, endcapEcalHitsToken_);
0098
0099
0100 if (thePhoton.r9() > minR9) {
0101
0102 double deltaE = scEnergyFunction_->getValue(*(thePhoton.superCluster()), 1);
0103 float e5x5 = thePhoton.e5x5();
0104 if (subdet == EcalBarrel)
0105 e5x5 = e5x5 * (1.0 + deltaE / thePhoton.superCluster()->rawEnergy());
0106 phoEcalEnergy = e5x5 + thePhoton.superCluster()->preshowerEnergy();
0107 } else {
0108 phoEcalEnergy = thePhoton.superCluster()->energy();
0109 }
0110
0111 thePhoton.setCorrectedEnergy(reco::Photon::ecal_standard, phoEcalEnergy, phoEcalEnergyError, false);
0112
0113
0114
0115 if (thePhoton.r9() > minR9) {
0116
0117 double deltaE = scEnergyFunction_->getValue(*(thePhoton.superCluster()), 1);
0118 float e5x5 = thePhoton.e5x5();
0119 if (subdet == EcalBarrel)
0120 e5x5 = e5x5 * (1.0 + deltaE / thePhoton.superCluster()->rawEnergy());
0121 phoEcalEnergy = e5x5 + thePhoton.superCluster()->preshowerEnergy();
0122
0123 phoEcalEnergy *= scCrackEnergyFunction_->getValue(*(thePhoton.superCluster()));
0124 phoEcalEnergyError = photonUncertaintyCalculator_->computePhotonEnergyUncertainty_highR9(
0125 thePhoton.superCluster()->eta(),
0126 thePhoton.superCluster()->phiWidth() / thePhoton.superCluster()->etaWidth(),
0127 phoEcalEnergy);
0128 } else {
0129
0130 phoEcalEnergy = photonEcalEnergyCorrFunction_->getValue(*(thePhoton.superCluster()), 1);
0131 phoEcalEnergy *= applyCrackCorrection(*(thePhoton.superCluster()), scCrackEnergyFunction_.get());
0132 phoEcalEnergyError = photonUncertaintyCalculator_->computePhotonEnergyUncertainty_lowR9(
0133 thePhoton.superCluster()->eta(),
0134 thePhoton.superCluster()->phiWidth() / thePhoton.superCluster()->etaWidth(),
0135 phoEcalEnergy);
0136 }
0137
0138
0139 thePhoton.setCorrectedEnergy(reco::Photon::ecal_photons, phoEcalEnergy, phoEcalEnergyError, false);
0140
0141
0142
0143 if ((weightsfromDB_ && !gedRegression_) || (!weightsfromDB_ && !(w_file_ == "none"))) {
0144 std::pair<double, double> cor =
0145 regressionCorrector_->CorrectedEnergyWithError(thePhoton, vtxcol, lazyTools, *theCaloGeom_);
0146 phoRegr1Energy = cor.first;
0147 phoRegr1EnergyError = cor.second;
0148
0149 thePhoton.setCorrectedEnergy(reco::Photon::regression1, phoRegr1Energy, phoRegr1EnergyError, false);
0150 }
0151
0152 if (gedRegression_) {
0153 gedRegression_->modifyObject(thePhoton);
0154
0155 thePhoton.setCorrectedEnergy(reco::Photon::regression1,
0156 thePhoton.getCorrectedEnergy(reco::Photon::regression2),
0157 thePhoton.getCorrectedEnergyError(reco::Photon::regression2),
0158 false);
0159 }
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172 }
0173
0174 double PhotonEnergyCorrector::applyCrackCorrection(const reco::SuperCluster& cl,
0175 EcalClusterFunctionBaseClass* crackCorrectionFunction) {
0176 double crackcor = 1.;
0177
0178 for (reco::CaloCluster_iterator cIt = cl.clustersBegin(); cIt != cl.clustersEnd(); ++cIt) {
0179 const reco::CaloClusterPtr cc = *cIt;
0180 crackcor *= ((cl.rawEnergy() + cc->energy() * (crackCorrectionFunction->getValue(*cc) - 1.)) / cl.rawEnergy());
0181 }
0182
0183 return crackcor;
0184 }