EgammaTools/src/PhotonEnergyCalibrator.cc

0001 #include "RecoEgamma/EgammaTools/interface/PhotonEnergyCalibrator.h"
0002
0003 #include "FWCore/ServiceRegistry/interface/Service.h"
0004 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007 #include <CLHEP/Random/RandGaussQ.h>
0008
0009 const EnergyScaleCorrection::ScaleCorrection PhotonEnergyCalibrator::defaultScaleCorr_;
0010 const EnergyScaleCorrection::SmearCorrection PhotonEnergyCalibrator::defaultSmearCorr_;
0011
0012 PhotonEnergyCalibrator::PhotonEnergyCalibrator(const std::string& correctionFile)
0013     : correctionRetriever_(correctionFile), rng_(nullptr), minEt_(1.0) {}
0014
0015 void PhotonEnergyCalibrator::initPrivateRng(TRandom* rnd) { rng_ = rnd; }
0016
0017 std::array<float, EGEnergySysIndex::kNrSysErrs> PhotonEnergyCalibrator::calibrate(
0018     reco::Photon& photon,
0019     const unsigned int runNumber,
0020     const EcalRecHitCollection* recHits,
0021     edm::StreamID const& id,
0022     const PhotonEnergyCalibrator::EventType eventType) const {
0023   return calibrate(photon, runNumber, recHits, gauss(id), eventType);
0024 }
0025
0026 std::array<float, EGEnergySysIndex::kNrSysErrs> PhotonEnergyCalibrator::calibrate(
0027     reco::Photon& photon,
0028     const unsigned int runNumber,
0029     const EcalRecHitCollection* recHits,
0030     const float smearNrSigma,
0031     const PhotonEnergyCalibrator::EventType eventType) const {
0032   const float scEtaAbs = std::abs(photon.superCluster()->eta());
0033   const float et = photon.getCorrectedEnergy(reco::Photon::P4type::regression2) / cosh(scEtaAbs);
0034
0035   if (et < minEt_ || edm::isNotFinite(et)) {
0036     std::array<float, EGEnergySysIndex::kNrSysErrs> retVal;
0037     retVal.fill(photon.getCorrectedEnergy(reco::Photon::P4type::regression2));
0038     retVal[EGEnergySysIndex::kScaleValue] = 1.0;
0039     retVal[EGEnergySysIndex::kSmearValue] = 0.0;
0040     retVal[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
0041     retVal[EGEnergySysIndex::kEcalErrPreCorr] = photon.getCorrectedEnergyError(reco::Photon::P4type::regression2);
0042     retVal[EGEnergySysIndex::kEcalErrPostCorr] = photon.getCorrectedEnergyError(reco::Photon::P4type::regression2);
0043     retVal[EGEnergySysIndex::kEcalTrkPreCorr] = 0.;
0044     retVal[EGEnergySysIndex::kEcalTrkErrPreCorr] = 0.;
0045     retVal[EGEnergySysIndex::kEcalTrkPostCorr] = 0.;
0046     retVal[EGEnergySysIndex::kEcalTrkErrPostCorr] = 0.;
0047
0048     return retVal;
0049   }
0050
0051   const DetId seedDetId = photon.superCluster()->seed()->seed();
0052   EcalRecHitCollection::const_iterator seedRecHit = recHits->find(seedDetId);
0053   unsigned int gainSeedSC = 12;
0054   if (seedRecHit != recHits->end()) {
0055     if (seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain6))
0056       gainSeedSC = 6;
0057     if (seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain1))
0058       gainSeedSC = 1;
0059   }
0060
0061   const EnergyScaleCorrection::ScaleCorrection* scaleCorr =
0062       correctionRetriever_.getScaleCorr(runNumber, et, scEtaAbs, photon.full5x5_r9(), gainSeedSC);
0063   const EnergyScaleCorrection::SmearCorrection* smearCorr =
0064       correctionRetriever_.getSmearCorr(runNumber, et, scEtaAbs, photon.full5x5_r9(), gainSeedSC);
0065   if (scaleCorr == nullptr)
0066     scaleCorr = &defaultScaleCorr_;
0067   if (smearCorr == nullptr)
0068     smearCorr = &defaultSmearCorr_;
0069
0070   std::array<float, EGEnergySysIndex::kNrSysErrs> uncertainties{};
0071
0072   uncertainties[EGEnergySysIndex::kScaleValue] = scaleCorr->scale();
0073   uncertainties[EGEnergySysIndex::kSmearValue] =
0074       smearCorr->sigma(et);  //even though we use scale = 1.0, we still store the value returned for MC
0075   uncertainties[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
0076   //MC central values are not scaled (scale = 1.0), data is not smeared (smearNrSigma = 0)
0077   //smearing still has a second order effect on data as it enters the E/p combination as an
0078   //extra uncertainty on the calo energy
0079   //MC gets all the scale systematics
0080   if (eventType == EventType::DATA) {
0081     setEnergyAndSystVarations(scaleCorr->scale(), 0., et, *scaleCorr, *smearCorr, photon, uncertainties);
0082   } else if (eventType == EventType::MC) {
0083     setEnergyAndSystVarations(1.0, smearNrSigma, et, *scaleCorr, *smearCorr, photon, uncertainties);
0084   }
0085
0086   return uncertainties;
0087 }
0088
0089 void PhotonEnergyCalibrator::setEnergyAndSystVarations(
0090     const float scale,
0091     const float smearNrSigma,
0092     const float et,
0093     const EnergyScaleCorrection::ScaleCorrection& scaleCorr,
0094     const EnergyScaleCorrection::SmearCorrection& smearCorr,
0095     reco::Photon& photon,
0096     std::array<float, EGEnergySysIndex::kNrSysErrs>& energyData) const {
0097   const float smear = smearCorr.sigma(et);
0098   const float smearRhoUp = smearCorr.sigma(et, 1, 0);
0099   const float smearRhoDn = smearCorr.sigma(et, -1, 0);
0100   const float smearPhiUp = smearCorr.sigma(et, 0, 1);
0101   const float smearPhiDn = smearCorr.sigma(et, 0, -1);
0102
0103   const float corr = scale + smear * smearNrSigma;
0104   const float corrRhoUp = scale + smearRhoUp * smearNrSigma;
0105   const float corrRhoDn = scale + smearRhoDn * smearNrSigma;
0106   const float corrPhiUp = scale + smearPhiUp * smearNrSigma;
0107   const float corrPhiDn = scale + smearPhiDn * smearNrSigma;
0108   const float corrUp = corrRhoUp;
0109   const float corrDn = corrRhoDn;
0110
0111   const double oldEcalEnergy = photon.getCorrectedEnergy(reco::Photon::P4type::regression2);
0112   const double oldEcalEnergyError = photon.getCorrectedEnergyError(reco::Photon::P4type::regression2);
0113
0114   energyData[EGEnergySysIndex::kEcalPreCorr] = oldEcalEnergy;
0115   energyData[EGEnergySysIndex::kEcalErrPreCorr] = oldEcalEnergyError;
0116
0117   const double newEcalEnergy = oldEcalEnergy * corr;
0118   const double newEcalEnergyError = std::hypot(oldEcalEnergyError * corr, smear * newEcalEnergy);
0119   photon.setCorrectedEnergy(reco::Photon::P4type::regression2, newEcalEnergy, newEcalEnergyError, true);
0120
0121   energyData[EGEnergySysIndex::kScaleStatUp] = oldEcalEnergy * (corr + scaleCorr.scaleErrStat());
0122   energyData[EGEnergySysIndex::kScaleStatDown] = oldEcalEnergy * (corr - scaleCorr.scaleErrStat());
0123   energyData[EGEnergySysIndex::kScaleSystUp] = oldEcalEnergy * (corr + scaleCorr.scaleErrSyst());
0124   energyData[EGEnergySysIndex::kScaleSystDown] = oldEcalEnergy * (corr - scaleCorr.scaleErrSyst());
0125   energyData[EGEnergySysIndex::kScaleGainUp] = oldEcalEnergy * (corr + scaleCorr.scaleErrGain());
0126   energyData[EGEnergySysIndex::kScaleGainDown] = oldEcalEnergy * (corr - scaleCorr.scaleErrGain());
0127   energyData[EGEnergySysIndex::kSmearRhoUp] = oldEcalEnergy * corrRhoUp;
0128   energyData[EGEnergySysIndex::kSmearRhoDown] = oldEcalEnergy * corrRhoDn;
0129   energyData[EGEnergySysIndex::kSmearPhiUp] = oldEcalEnergy * corrPhiUp;
0130   energyData[EGEnergySysIndex::kSmearPhiDown] = oldEcalEnergy * corrPhiDn;
0131
0132   // The total variation
0133   energyData[EGEnergySysIndex::kScaleUp] =
0134       oldEcalEnergy * (corr + scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain));
0135   energyData[EGEnergySysIndex::kScaleDown] =
0136       oldEcalEnergy * (corr - scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain));
0137   energyData[EGEnergySysIndex::kSmearUp] = oldEcalEnergy * corrUp;
0138   energyData[EGEnergySysIndex::kSmearDown] = oldEcalEnergy * corrDn;
0139
0140   energyData[EGEnergySysIndex::kEcalPostCorr] = photon.getCorrectedEnergy(reco::Photon::P4type::regression2);
0141   energyData[EGEnergySysIndex::kEcalErrPostCorr] = photon.getCorrectedEnergyError(reco::Photon::P4type::regression2);
0142 }
0143
0144 double PhotonEnergyCalibrator::gauss(edm::StreamID const& id) const {
0145   if (rng_) {
0146     return rng_->Gaus();
0147   } else {
0148     edm::Service<edm::RandomNumberGenerator> rng;
0149     if (!rng.isAvailable()) {
0150       throw cms::Exception("Configuration")
0151           << "XXXXXXX requires the RandomNumberGeneratorService\n"
0152              "which is not present in the configuration file.  You must add the service\n"
0153              "in the configuration file or remove the modules that require it.";
0154     }
0155     CLHEP::RandGaussQ gaussDistribution(rng->getEngine(id), 0.0, 1.0);
0156     return gaussDistribution.fire();
0157   }
0158 }