ElectronTools/src/PhotonEnergyCalibratorRun2.cc

0001 #include "EgammaAnalysis/ElectronTools/interface/PhotonEnergyCalibratorRun2.h"
0002 #include <CLHEP/Random/RandGaussQ.h>
0003 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0004 #include "FWCore/ServiceRegistry/interface/Service.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 PhotonEnergyCalibratorRun2::PhotonEnergyCalibratorRun2(bool isMC, bool synchronization, std::string correctionFile)
0008     : isMC_(isMC),
0009       synchronization_(synchronization),
0010       rng_(nullptr),
0011       _correctionRetriever(correctionFile)  // here is opening the files and reading thecorrections
0012 {
0013   if (isMC_) {
0014     _correctionRetriever.doScale = false;
0015     _correctionRetriever.doSmearings = true;
0016   } else {
0017     _correctionRetriever.doScale = true;
0018     _correctionRetriever.doSmearings = false;
0019   }
0020 }
0021
0022 PhotonEnergyCalibratorRun2::~PhotonEnergyCalibratorRun2() {}
0023
0024 void PhotonEnergyCalibratorRun2::initPrivateRng(TRandom *rnd) { rng_ = rnd; }
0025
0026 void PhotonEnergyCalibratorRun2::calibrate(reco::Photon &photon,
0027                                            unsigned int runNumber,
0028                                            edm::StreamID const &id) const {
0029   SimplePhoton simple(photon, runNumber, isMC_);
0030   calibrate(simple, id);
0031   simple.writeTo(photon);
0032 }
0033
0034 void PhotonEnergyCalibratorRun2::calibrate(SimplePhoton &photon, edm::StreamID const &id) const {
0035   assert(isMC_ == photon.isMC());
0036   float smear = 0.0, scale = 1.0;
0037   float aeta = std::abs(photon.getEta());  //, r9 = photon.getR9();
0038   float et = photon.getNewEnergy() / cosh(aeta);
0039
0040   scale = _correctionRetriever.ScaleCorrection(photon.getRunNumber(), photon.isEB(), photon.getR9(), aeta, et);
0041   smear = _correctionRetriever.getSmearingSigma(photon.getRunNumber(), photon.isEB(), photon.getR9(), aeta, et, 0., 0.);
0042
0043   double newEcalEnergy, newEcalEnergyError;
0044   if (isMC_) {
0045     double corr = 1.0 + smear * gauss(id);
0046     newEcalEnergy = photon.getNewEnergy() * corr;
0047     newEcalEnergyError = std::hypot(photon.getNewEnergyError() * corr, smear * newEcalEnergy);
0048   } else {
0049     newEcalEnergy = photon.getNewEnergy() * scale;
0050     newEcalEnergyError = std::hypot(photon.getNewEnergyError() * scale, smear * newEcalEnergy);
0051   }
0052   photon.setNewEnergy(newEcalEnergy);
0053   photon.setNewEnergyError(newEcalEnergyError);
0054 }
0055
0056 double PhotonEnergyCalibratorRun2::gauss(edm::StreamID const &id) const {
0057   if (synchronization_)
0058     return 1.0;
0059   if (rng_) {
0060     return rng_->Gaus();
0061   } else {
0062     edm::Service<edm::RandomNumberGenerator> rng;
0063     if (!rng.isAvailable()) {
0064       throw cms::Exception("Configuration")
0065           << "XXXXXXX requires the RandomNumberGeneratorService\n"
0066              "which is not present in the configuration file.  You must add the service\n"
0067              "in the configuration file or remove the modules that require it.";
0068     }
0069     CLHEP::RandGaussQ gaussDistribution(rng->getEngine(id), 0.0, 1.0);
0070     return gaussDistribution.fire();
0071   }
0072 }