EgammaTools/plugins/CalibratedPhotonProducers.cc

0001 //author: Alan Smithee
0002 //description:
0003 //  this class allows the residual scale and smearing to be applied to photon
0004 //  it will write out all the calibration info in the event, such as scale correction value,
0005 //  smearing correction value, random nr used, energy post calibration, energy pre calibration
0006 //  can optionally write out a new collection of photon with the energy corrected by default
0007 //  a port of EgammaAnalysis/ElectronTools/CalibratedPhotonProducerRun2
0008
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "FWCore/Framework/interface/stream/EDProducer.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/Utilities/interface/EDGetToken.h"
0018
0019 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0020 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0021 #include "DataFormats/PatCandidates/interface/Photon.h"
0022 #include "RecoEgamma/EgammaTools/interface/PhotonEnergyCalibrator.h"
0023 #include "RecoEgamma/EgammaTools/interface/EGEnergySysIndex.h"
0024 #include "RecoEgamma/EgammaTools/interface/EgammaRandomSeeds.h"
0025
0026 #include "TRandom2.h"
0027
0028 #include <memory>
0029
0030 #include <random>
0031 #include <vector>
0032
0033 template <typename T>
0034 class CalibratedPhotonProducerT : public edm::stream::EDProducer<> {
0035 public:
0036   explicit CalibratedPhotonProducerT(const edm::ParameterSet&);
0037   ~CalibratedPhotonProducerT() override {}
0038   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0039   void produce(edm::Event&, const edm::EventSetup&) override;
0040
0041 private:
0042   void setSemiDetRandomSeed(const edm::Event& iEvent, const T& obj, size_t nrObjs, size_t objNr);
0043
0044   edm::EDGetTokenT<edm::View<T>> photonToken_;
0045   PhotonEnergyCalibrator energyCorrector_;
0046   std::unique_ptr<TRandom> semiDeterministicRng_;
0047   edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEBToken_;
0048   edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEEToken_;
0049   bool produceCalibratedObjs_;
0050
0051   static const std::vector<int> valMapsToStore_;
0052 };
0053
0054 template <typename T>
0055 const std::vector<int> CalibratedPhotonProducerT<T>::valMapsToStore_ = {
0056     EGEnergySysIndex::kScaleStatUp,    EGEnergySysIndex::kScaleStatDown, EGEnergySysIndex::kScaleSystUp,
0057     EGEnergySysIndex::kScaleSystDown,  EGEnergySysIndex::kScaleGainUp,   EGEnergySysIndex::kScaleGainDown,
0058     EGEnergySysIndex::kSmearRhoUp,     EGEnergySysIndex::kSmearRhoDown,  EGEnergySysIndex::kSmearPhiUp,
0059     EGEnergySysIndex::kSmearPhiDown,   EGEnergySysIndex::kScaleUp,       EGEnergySysIndex::kScaleDown,
0060     EGEnergySysIndex::kSmearUp,        EGEnergySysIndex::kSmearDown,     EGEnergySysIndex::kScaleValue,
0061     EGEnergySysIndex::kSmearValue,     EGEnergySysIndex::kSmearNrSigma,  EGEnergySysIndex::kEcalPreCorr,
0062     EGEnergySysIndex::kEcalErrPreCorr, EGEnergySysIndex::kEcalPostCorr,  EGEnergySysIndex::kEcalErrPostCorr};
0063
0064 namespace {
0065   template <typename HandleType, typename ValType>
0066   void fillAndStoreValueMap(edm::Event& iEvent,
0067                             HandleType objHandle,
0068                             const std::vector<ValType>& vals,
0069                             const std::string& name) {
0070     auto valMap = std::make_unique<edm::ValueMap<ValType>>();
0071     typename edm::ValueMap<ValType>::Filler filler(*valMap);
0072     filler.insert(objHandle, vals.begin(), vals.end());
0073     filler.fill();
0074     iEvent.put(std::move(valMap), name);
0075   }
0076 }  // namespace
0077
0078 template <typename T>
0079 CalibratedPhotonProducerT<T>::CalibratedPhotonProducerT(const edm::ParameterSet& conf)
0080     : photonToken_(consumes(conf.getParameter<edm::InputTag>("src"))),
0081       energyCorrector_(conf.getParameter<std::string>("correctionFile")),
0082       recHitCollectionEBToken_(consumes(conf.getParameter<edm::InputTag>("recHitCollectionEB"))),
0083       recHitCollectionEEToken_(consumes(conf.getParameter<edm::InputTag>("recHitCollectionEE"))),
0084       produceCalibratedObjs_(conf.getParameter<bool>("produceCalibratedObjs")) {
0085   energyCorrector_.setMinEt(conf.getParameter<double>("minEtToCalibrate"));
0086
0087   if (conf.getParameter<bool>("semiDeterministic")) {
0088     semiDeterministicRng_ = std::make_unique<TRandom2>();
0089     energyCorrector_.initPrivateRng(semiDeterministicRng_.get());
0090   }
0091
0092   if (produceCalibratedObjs_)
0093     produces<std::vector<T>>();
0094
0095   for (const auto& toStore : valMapsToStore_) {
0096     produces<edm::ValueMap<float>>(EGEnergySysIndex::name(toStore));
0097   }
0098 }
0099
0100 template <typename T>
0101 void CalibratedPhotonProducerT<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0102   edm::ParameterSetDescription desc;
0103   desc.add<edm::InputTag>("src", edm::InputTag("gedGsfElectrons"));
0104   desc.add<edm::InputTag>("recHitCollectionEB", edm::InputTag("reducedEcalRecHitsEB"));
0105   desc.add<edm::InputTag>("recHitCollectionEE", edm::InputTag("reducedEcalRecHitsEE"));
0106   desc.add<std::string>("correctionFile", std::string());
0107   desc.add<double>("minEtToCalibrate", 5.0);
0108   desc.add<bool>("produceCalibratedObjs", true);
0109   desc.add<bool>("semiDeterministic", true);
0110   std::vector<std::string> valMapsProduced;
0111   valMapsProduced.reserve(valMapsToStore_.size());
0112   for (auto varToStore : valMapsToStore_)
0113     valMapsProduced.push_back(EGEnergySysIndex::name(varToStore));
0114   desc.add<std::vector<std::string>>("valueMapsStored", valMapsProduced)
0115       ->setComment(
0116           "provides to python configs the list of valuemaps stored, can not be overriden in the python config");
0117   descriptions.addWithDefaultLabel(desc);
0118 }
0119
0120 template <typename T>
0121 void CalibratedPhotonProducerT<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0122   auto inHandle = iEvent.getHandle(photonToken_);
0123
0124   auto recHitCollectionEBHandle = iEvent.getHandle(recHitCollectionEBToken_);
0125   auto recHitCollectionEEHandle = iEvent.getHandle(recHitCollectionEEToken_);
0126
0127   std::unique_ptr<std::vector<T>> out = std::make_unique<std::vector<T>>();
0128
0129   size_t nrObj = inHandle->size();
0130   std::array<std::vector<float>, EGEnergySysIndex::kNrSysErrs> results;
0131   for (auto& res : results)
0132     res.reserve(nrObj);
0133
0134   const PhotonEnergyCalibrator::EventType evtType =
0135       iEvent.isRealData() ? PhotonEnergyCalibrator::EventType::DATA : PhotonEnergyCalibrator::EventType::MC;
0136
0137   for (const auto& pho : *inHandle) {
0138     out->emplace_back(pho);
0139
0140     if (semiDeterministicRng_)
0141       setSemiDetRandomSeed(iEvent, pho, nrObj, out->size());
0142
0143     const EcalRecHitCollection* recHits =
0144         (pho.isEB()) ? recHitCollectionEBHandle.product() : recHitCollectionEEHandle.product();
0145     std::array<float, EGEnergySysIndex::kNrSysErrs> uncertainties =
0146         energyCorrector_.calibrate(out->back(), iEvent.id().run(), recHits, iEvent.streamID(), evtType);
0147
0148     for (size_t index = 0; index < EGEnergySysIndex::kNrSysErrs; index++) {
0149       results[index].push_back(uncertainties[index]);
0150     }
0151   }
0152
0153   auto fillAndStore = [&](auto handle) {
0154     for (const auto& mapToStore : valMapsToStore_) {
0155       fillAndStoreValueMap(iEvent, handle, results[mapToStore], EGEnergySysIndex::name(mapToStore));
0156     }
0157   };
0158
0159   if (produceCalibratedObjs_) {
0160     fillAndStore(iEvent.put(std::move(out)));
0161   } else {
0162     fillAndStore(inHandle);
0163   }
0164 }
0165
0166 //needs to be synced to CalibratedElectronProducers, want the same seed for a given SC
0167 template <typename T>
0168 void CalibratedPhotonProducerT<T>::setSemiDetRandomSeed(const edm::Event& iEvent,
0169                                                         const T& obj,
0170                                                         size_t nrObjs,
0171                                                         size_t objNr) {
0172   if (obj.superCluster().isNonnull()) {
0173     semiDeterministicRng_->SetSeed(egamma::getRandomSeedFromSC(iEvent, obj.superCluster()));
0174   } else {
0175     semiDeterministicRng_->SetSeed(egamma::getRandomSeedFromObj(iEvent, obj, nrObjs, objNr));
0176   }
0177 }
0178
0179 using CalibratedPhotonProducer = CalibratedPhotonProducerT<reco::Photon>;
0180 using CalibratedPatPhotonProducer = CalibratedPhotonProducerT<pat::Photon>;
0181
0182 #include "FWCore/Framework/interface/MakerMacros.h"
0183
0184 DEFINE_FWK_MODULE(CalibratedPhotonProducer);
0185 DEFINE_FWK_MODULE(CalibratedPatPhotonProducer);