EgammaTools/src/ElectronEnergyCalibrator.cc

0001 #include "RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.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 ElectronEnergyCalibrator::defaultScaleCorr_;
0010 const EnergyScaleCorrection::SmearCorrection ElectronEnergyCalibrator::defaultSmearCorr_;
0011
0012 ElectronEnergyCalibrator::ElectronEnergyCalibrator(const EpCombinationTool& combinator,
0013                                                    const std::string& correctionFile)
0014     : correctionRetriever_(correctionFile), epCombinationTool_(&combinator), rng_(nullptr), minEt_(1.0) {}
0015
0016 void ElectronEnergyCalibrator::initPrivateRng(TRandom* rnd) { rng_ = rnd; }
0017
0018 std::array<float, EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::calibrate(
0019     reco::GsfElectron& ele,
0020     const unsigned int runNumber,
0021     const EcalRecHitCollection* recHits,
0022     edm::StreamID const& id,
0023     const ElectronEnergyCalibrator::EventType eventType) const {
0024   return calibrate(ele, runNumber, recHits, gauss(id), eventType);
0025 }
0026
0027 std::array<float, EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::calibrate(
0028     reco::GsfElectron& ele,
0029     unsigned int runNumber,
0030     const EcalRecHitCollection* recHits,
0031     const float smearNrSigma,
0032     const ElectronEnergyCalibrator::EventType eventType) const {
0033   const float scEtaAbs = std::abs(ele.superCluster()->eta());
0034   const float et = ele.ecalEnergy() / cosh(scEtaAbs);
0035
0036   if (et < minEt_ || edm::isNotFinite(et)) {
0037     std::array<float, EGEnergySysIndex::kNrSysErrs> retVal;
0038     retVal.fill(ele.energy());
0039     retVal[EGEnergySysIndex::kScaleValue] = 1.0;
0040     retVal[EGEnergySysIndex::kSmearValue] = 0.0;
0041     retVal[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
0042     retVal[EGEnergySysIndex::kEcalPreCorr] = ele.ecalEnergy();
0043     retVal[EGEnergySysIndex::kEcalErrPreCorr] = ele.ecalEnergyError();
0044     retVal[EGEnergySysIndex::kEcalPostCorr] = ele.ecalEnergy();
0045     retVal[EGEnergySysIndex::kEcalErrPostCorr] = ele.ecalEnergyError();
0046     retVal[EGEnergySysIndex::kEcalTrkPreCorr] = ele.energy();
0047     retVal[EGEnergySysIndex::kEcalTrkErrPreCorr] = ele.corrections().combinedP4Error;
0048     retVal[EGEnergySysIndex::kEcalTrkPostCorr] = ele.energy();
0049     retVal[EGEnergySysIndex::kEcalTrkErrPostCorr] = ele.corrections().combinedP4Error;
0050     return retVal;
0051   }
0052
0053   const DetId seedDetId = ele.superCluster()->seed()->seed();
0054   EcalRecHitCollection::const_iterator seedRecHit = recHits->find(seedDetId);
0055   unsigned int gainSeedSC = 12;
0056   if (seedRecHit != recHits->end()) {
0057     if (seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain6))
0058       gainSeedSC = 6;
0059     if (seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain1))
0060       gainSeedSC = 1;
0061   }
0062
0063   const EnergyScaleCorrection::ScaleCorrection* scaleCorr =
0064       correctionRetriever_.getScaleCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
0065   const EnergyScaleCorrection::SmearCorrection* smearCorr =
0066       correctionRetriever_.getSmearCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
0067   if (scaleCorr == nullptr)
0068     scaleCorr = &defaultScaleCorr_;
0069   if (smearCorr == nullptr)
0070     smearCorr = &defaultSmearCorr_;
0071
0072   std::array<float, EGEnergySysIndex::kNrSysErrs> uncertainties{};
0073
0074   uncertainties[EGEnergySysIndex::kScaleValue] = scaleCorr->scale();
0075   uncertainties[EGEnergySysIndex::kSmearValue] =
0076       smearCorr->sigma(et);  //even though we use scale = 1.0, we still store the value returned for MC
0077   uncertainties[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
0078   //MC central values are not scaled (scale = 1.0), data is not smeared (smearNrSigma = 0)
0079   //the smearing (or resolution extra parameter as it might better be called)
0080   //still has a second order effect on data as it enters the E/p combination as an adjustment
0081   //to the estimate of the resolution contained in caloEnergyError
0082   //MC gets all the scale systematics
0083   if (eventType == EventType::DATA) {
0084     setEnergyAndSystVarations(scaleCorr->scale(), 0., et, *scaleCorr, *smearCorr, ele, uncertainties);
0085   } else if (eventType == EventType::MC) {
0086     setEnergyAndSystVarations(1.0, smearNrSigma, et, *scaleCorr, *smearCorr, ele, uncertainties);
0087   }
0088
0089   return uncertainties;
0090 }
0091
0092 void ElectronEnergyCalibrator::setEnergyAndSystVarations(
0093     const float scale,
0094     const float smearNrSigma,
0095     const float et,
0096     const EnergyScaleCorrection::ScaleCorrection& scaleCorr,
0097     const EnergyScaleCorrection::SmearCorrection& smearCorr,
0098     reco::GsfElectron& ele,
0099     std::array<float, EGEnergySysIndex::kNrSysErrs>& energyData) const {
0100   const float smear = smearCorr.sigma(et);
0101   const float smearRhoUp = smearCorr.sigma(et, 1, 0);
0102   const float smearRhoDn = smearCorr.sigma(et, -1, 0);
0103   const float smearPhiUp = smearCorr.sigma(et, 0, 1);
0104   const float smearPhiDn = smearCorr.sigma(et, 0, -1);
0105   const float smearUp = smearRhoUp;
0106   const float smearDn = smearRhoDn;
0107
0108   const float corr = scale + smear * smearNrSigma;
0109   const float corrRhoUp = scale + smearRhoUp * smearNrSigma;
0110   const float corrRhoDn = scale + smearRhoDn * smearNrSigma;
0111   const float corrPhiUp = scale + smearPhiUp * smearNrSigma;
0112   const float corrPhiDn = scale + smearPhiDn * smearNrSigma;
0113   const float corrUp = corrRhoUp;
0114   const float corrDn = corrRhoDn;
0115
0116   const float corrScaleStatUp = corr + scaleCorr.scaleErrStat();
0117   const float corrScaleStatDn = corr - scaleCorr.scaleErrStat();
0118   const float corrScaleSystUp = corr + scaleCorr.scaleErrSyst();
0119   const float corrScaleSystDn = corr - scaleCorr.scaleErrSyst();
0120   const float corrScaleGainUp = corr + scaleCorr.scaleErrGain();
0121   const float corrScaleGainDn = corr - scaleCorr.scaleErrGain();
0122   const float corrScaleUp = corr + scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
0123   const float corrScaleDn = corr - scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
0124
0125   const math::XYZTLorentzVector oldP4 = ele.p4();
0126   energyData[EGEnergySysIndex::kEcalTrkPreCorr] = ele.energy();
0127   energyData[EGEnergySysIndex::kEcalTrkErrPreCorr] = ele.corrections().combinedP4Error;
0128   energyData[EGEnergySysIndex::kEcalPreCorr] = ele.ecalEnergy();
0129   energyData[EGEnergySysIndex::kEcalErrPreCorr] = ele.ecalEnergyError();
0130
0131   energyData[EGEnergySysIndex::kScaleStatUp] = calCombinedMom(ele, corrScaleStatUp, smear).first;
0132   energyData[EGEnergySysIndex::kScaleStatDown] = calCombinedMom(ele, corrScaleStatDn, smear).first;
0133   energyData[EGEnergySysIndex::kScaleSystUp] = calCombinedMom(ele, corrScaleSystUp, smear).first;
0134   energyData[EGEnergySysIndex::kScaleSystDown] = calCombinedMom(ele, corrScaleSystDn, smear).first;
0135   energyData[EGEnergySysIndex::kScaleGainUp] = calCombinedMom(ele, corrScaleGainUp, smear).first;
0136   energyData[EGEnergySysIndex::kScaleGainDown] = calCombinedMom(ele, corrScaleGainDn, smear).first;
0137
0138   energyData[EGEnergySysIndex::kSmearRhoUp] = calCombinedMom(ele, corrRhoUp, smearRhoUp).first;
0139   energyData[EGEnergySysIndex::kSmearRhoDown] = calCombinedMom(ele, corrRhoDn, smearRhoDn).first;
0140   energyData[EGEnergySysIndex::kSmearPhiUp] = calCombinedMom(ele, corrPhiUp, smearPhiUp).first;
0141   energyData[EGEnergySysIndex::kSmearPhiDown] = calCombinedMom(ele, corrPhiDn, smearPhiDn).first;
0142
0143   energyData[EGEnergySysIndex::kScaleUp] = calCombinedMom(ele, corrScaleUp, smear).first;
0144   energyData[EGEnergySysIndex::kScaleDown] = calCombinedMom(ele, corrScaleDn, smear).first;
0145   energyData[EGEnergySysIndex::kSmearUp] = calCombinedMom(ele, corrUp, smearUp).first;
0146   energyData[EGEnergySysIndex::kSmearDown] = calCombinedMom(ele, corrDn, smearDn).first;
0147
0148   const std::pair<float, float> combinedMomentum = calCombinedMom(ele, corr, smear);
0149   setEcalEnergy(ele, corr, smear);
0150   const float energyCorr = combinedMomentum.first / oldP4.t();
0151
0152   const math::XYZTLorentzVector newP4(
0153       oldP4.x() * energyCorr, oldP4.y() * energyCorr, oldP4.z() * energyCorr, combinedMomentum.first);
0154
0155   ele.correctMomentum(newP4, ele.trackMomentumError(), combinedMomentum.second);
0156   energyData[EGEnergySysIndex::kEcalTrkPostCorr] = ele.energy();
0157   energyData[EGEnergySysIndex::kEcalTrkErrPostCorr] = ele.corrections().combinedP4Error;
0158
0159   energyData[EGEnergySysIndex::kEcalPostCorr] = ele.ecalEnergy();
0160   energyData[EGEnergySysIndex::kEcalErrPostCorr] = ele.ecalEnergyError();
0161 }
0162
0163 void ElectronEnergyCalibrator::setEcalEnergy(reco::GsfElectron& ele, const float scale, const float smear) const {
0164   const float oldEcalEnergy = ele.ecalEnergy();
0165   const float oldEcalEnergyErr = ele.ecalEnergyError();
0166   ele.setCorrectedEcalEnergy(oldEcalEnergy * scale);
0167   ele.setCorrectedEcalEnergyError(std::hypot(oldEcalEnergyErr * scale, oldEcalEnergy * smear * scale));
0168 }
0169
0170 std::pair<float, float> ElectronEnergyCalibrator::calCombinedMom(reco::GsfElectron& ele,
0171                                                                  const float scale,
0172                                                                  const float smear) const {
0173   const float oldEcalEnergy = ele.ecalEnergy();
0174   const float oldEcalEnergyErr = ele.ecalEnergyError();
0175
0176   const auto oldP4 = ele.p4();
0177   const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
0178   const float oldTrkMomErr = ele.trackMomentumError();
0179
0180   setEcalEnergy(ele, scale, smear);
0181   const auto& combinedMomentum = epCombinationTool_->combine(ele, oldEcalEnergyErr * scale);
0182
0183   ele.setCorrectedEcalEnergy(oldEcalEnergy);
0184   ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
0185   ele.correctMomentum(oldP4, oldTrkMomErr, oldP4Err);
0186
0187   return combinedMomentum;
0188 }
0189
0190 double ElectronEnergyCalibrator::gauss(edm::StreamID const& id) const {
0191   if (rng_) {
0192     return rng_->Gaus();
0193   } else {
0194     edm::Service<edm::RandomNumberGenerator> rng;
0195     if (!rng.isAvailable()) {
0196       throw cms::Exception("Configuration")
0197           << "XXXXXXX requires the RandomNumberGeneratorService\n"
0198              "which is not present in the configuration file.  You must add the service\n"
0199              "in the configuration file or remove the modules that require it.";
0200     }
0201     CLHEP::RandGaussQ gaussDistribution(rng->getEngine(id), 0.0, 1.0);
0202     return gaussDistribution.fire();
0203   }
0204 }