EgammaPhotonAlgos/src/PhotonEnergyCorrector.cc

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   // get the geometry from the event setup:
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   //  candidateP4type_ = config.getParameter<std::string>("candidateP4type") ;
0026
0027   // function to extract f(eta) correction
0028   std::string superClusterFunctionName = config.getParameter<std::string>("superClusterEnergyCorrFunction");
0029   scEnergyFunction_ = EcalClusterFunctionFactory::get()->create(superClusterFunctionName, config, iC);
0030
0031   // function to extract corrections to cracks
0032   std::string superClusterCrackFunctionName = config.getParameter<std::string>("superClusterCrackEnergyCorrFunction");
0033   scCrackEnergyFunction_ = EcalClusterFunctionFactory::get()->create(superClusterCrackFunctionName, config, iC);
0034
0035   // function to extract the error on the sc ecal correction
0036   std::string superClusterErrorFunctionName = config.getParameter<std::string>("superClusterEnergyErrorFunction");
0037   scEnergyErrorFunction_ = EcalClusterFunctionFactory::get()->create(superClusterErrorFunctionName, config, iC);
0038
0039   // function  to extract the error on the photon ecal correction
0040   std::string photonEnergyFunctionName = config.getParameter<std::string>("photonEcalEnergyCorrFunction");
0041   photonEcalEnergyCorrFunction_ = EcalClusterFunctionFactory::get()->create(photonEnergyFunctionName, config, iC);
0042   //ingredient for photon uncertainty
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   // ingredient for energy regression
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   ////////////// Here default Ecal corrections based on electrons  ////////////////////////
0100   if (thePhoton.r9() > minR9) {
0101     // f(eta) correction to e5x5
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   // store the value in the Photon.h
0111   thePhoton.setCorrectedEnergy(reco::Photon::ecal_standard, phoEcalEnergy, phoEcalEnergyError, false);
0112
0113   ////////////// Here Ecal corrections specific for photons ////////////////////////
0114
0115   if (thePhoton.r9() > minR9) {
0116     // f(eta) correction to e5x5
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     // add correction for cracks
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     // correction for low r9
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   // store the value in the Photon.h
0139   thePhoton.setCorrectedEnergy(reco::Photon::ecal_photons, phoEcalEnergy, phoEcalEnergyError, false);
0140
0141   //////////  Energy  Regression //////////////////////
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     // store the value in the Photon.h
0149     thePhoton.setCorrectedEnergy(reco::Photon::regression1, phoRegr1Energy, phoRegr1EnergyError, false);
0150   }
0151
0152   if (gedRegression_) {
0153     gedRegression_->modifyObject(thePhoton);  // uses regression2 slot
0154     // force regresions1 and 2 to be the same (no reason to be different)
0155     thePhoton.setCorrectedEnergy(reco::Photon::regression1,
0156                                  thePhoton.getCorrectedEnergy(reco::Photon::regression2),
0157                                  thePhoton.getCorrectedEnergyError(reco::Photon::regression2),
0158                                  false);
0159   }
0160
0161   /*
0162   std::cout << " ------------------------- " << std::endl;
0163   std::cout << " Corrector " << std::endl;
0164   std::cout << " P4 Type " << thePhoton.getCandidateP4type() << " candidate p4 " << thePhoton.p4() << std::endl;
0165   std::cout << " photon ecalEnergy " << thePhoton.getCorrectedEnergy(reco::Photon::ecal_photons) << " error " << thePhoton.getCorrectedEnergyError(reco::Photon::ecal_photons) << std::endl;
0166   std::cout << " ecal p4 from accessor " << thePhoton.p4(reco::Photon::ecal_photons) <<  std::endl;
0167   std::cout << " ------------------------- " << std::endl;
0168   std::cout << " reg1 energy " << thePhoton.getCorrectedEnergy(reco::Photon::regression1)  << " error " <<  thePhoton.getCorrectedEnergyError(reco::Photon::regression1) << std::endl;
0169   std::cout << " New p4 from regression " <<  thePhoton.p4(reco::Photon::regression1)    << std::endl;
0170   std::cout << " ------------------------- " << std::endl;
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   }  // loop on BCs
0182
0183   return crackcor;
0184 }