CalibratedPhotonProducers.cc

Go to the documentation of this file.

//author: Alan Smithee
//description:
//  this class allows the residual scale and smearing to be applied to photon
//  it will write out all the calibration info in the event, such as scale correction value,
//  smearing correction value, random nr used, energy post calibration, energy pre calibration
//  can optionally write out a new collection of photon with the energy corrected by default
//  a port of EgammaAnalysis/ElectronTools/CalibratedPhotonProducerRun2

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "RecoEgamma/EgammaTools/interface/PhotonEnergyCalibrator.h"
#include "RecoEgamma/EgammaTools/interface/EGEnergySysIndex.h"
#include "RecoEgamma/EgammaTools/interface/EgammaRandomSeeds.h"

#include "TRandom2.h"

#include <vector>
#include <random>

template <typename T>
class CalibratedPhotonProducerT : public edm::stream::EDProducer<> {
public:
  explicit CalibratedPhotonProducerT(const edm::ParameterSet&);
  ~CalibratedPhotonProducerT() override{};
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
  void produce(edm::Event&, const edm::EventSetup&) override;

private:
  void setSemiDetRandomSeed(const edm::Event& iEvent, const T& obj, size_t nrObjs, size_t objNr);

  edm::EDGetTokenT<edm::View<T>> photonToken_;
  PhotonEnergyCalibrator energyCorrector_;
  std::unique_ptr<TRandom> semiDeterministicRng_;
  edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEBToken_;
  edm::EDGetTokenT<EcalRecHitCollection> recHitCollectionEEToken_;
  bool produceCalibratedObjs_;

  static const std::vector<int> valMapsToStore_;
};

template <typename T>
const std::vector<int> CalibratedPhotonProducerT<T>::valMapsToStore_ = {
    EGEnergySysIndex::kScaleStatUp,    EGEnergySysIndex::kScaleStatDown, EGEnergySysIndex::kScaleSystUp,
    EGEnergySysIndex::kScaleSystDown,  EGEnergySysIndex::kScaleGainUp,   EGEnergySysIndex::kScaleGainDown,
    EGEnergySysIndex::kSmearRhoUp,     EGEnergySysIndex::kSmearRhoDown,  EGEnergySysIndex::kSmearPhiUp,
    EGEnergySysIndex::kSmearPhiDown,   EGEnergySysIndex::kScaleUp,       EGEnergySysIndex::kScaleDown,
    EGEnergySysIndex::kSmearUp,        EGEnergySysIndex::kSmearDown,     EGEnergySysIndex::kScaleValue,
    EGEnergySysIndex::kSmearValue,     EGEnergySysIndex::kSmearNrSigma,  EGEnergySysIndex::kEcalPreCorr,
    EGEnergySysIndex::kEcalErrPreCorr, EGEnergySysIndex::kEcalPostCorr,  EGEnergySysIndex::kEcalErrPostCorr};

namespace {
  template <typename HandleType, typename ValType>
  void fillAndStoreValueMap(edm::Event& iEvent,
                            HandleType objHandle,
                            const std::vector<ValType>& vals,
                            const std::string& name) {
    auto valMap = std::make_unique<edm::ValueMap<ValType>>();
    typename edm::ValueMap<ValType>::Filler filler(*valMap);
    filler.insert(objHandle, vals.begin(), vals.end());
    filler.fill();
    iEvent.put(std::move(valMap), name);
  }
}  // namespace

template <typename T>
CalibratedPhotonProducerT<T>::CalibratedPhotonProducerT(const edm::ParameterSet& conf)
    : photonToken_(consumes<edm::View<T>>(conf.getParameter<edm::InputTag>("src"))),
      energyCorrector_(conf.getParameter<std::string>("correctionFile")),
      recHitCollectionEBToken_(consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("recHitCollectionEB"))),
      recHitCollectionEEToken_(consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("recHitCollectionEE"))),
      produceCalibratedObjs_(conf.getParameter<bool>("produceCalibratedObjs")) {
  energyCorrector_.setMinEt(conf.getParameter<double>("minEtToCalibrate"));

  if (conf.getParameter<bool>("semiDeterministic")) {
    semiDeterministicRng_.reset(new TRandom2());
    energyCorrector_.initPrivateRng(semiDeterministicRng_.get());
  }

  if (produceCalibratedObjs_)
    produces<std::vector<T>>();

  for (const auto& toStore : valMapsToStore_) {
    produces<edm::ValueMap<float>>(EGEnergySysIndex::name(toStore));
  }
}

template <typename T>
void CalibratedPhotonProducerT<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("src", edm::InputTag("gedGsfElectrons"));
  desc.add<edm::InputTag>("recHitCollectionEB", edm::InputTag("reducedEcalRecHitsEB"));
  desc.add<edm::InputTag>("recHitCollectionEE", edm::InputTag("reducedEcalRecHitsEE"));
  desc.add<std::string>("correctionFile", std::string());
  desc.add<double>("minEtToCalibrate", 5.0);
  desc.add<bool>("produceCalibratedObjs", true);
  desc.add<bool>("semiDeterministic", true);
  std::vector<std::string> valMapsProduced;
  for (auto varToStore : valMapsToStore_)
    valMapsProduced.push_back(EGEnergySysIndex::name(varToStore));
  desc.add<std::vector<std::string>>("valueMapsStored", valMapsProduced)
      ->setComment(
          "provides to python configs the list of valuemaps stored, can not be overriden in the python config");
  descriptions.addWithDefaultLabel(desc);
}

template <typename T>
void CalibratedPhotonProducerT<T>::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::Handle<edm::View<T>> inHandle;
  iEvent.getByToken(photonToken_, inHandle);

  edm::Handle<EcalRecHitCollection> recHitCollectionEBHandle;
  edm::Handle<EcalRecHitCollection> recHitCollectionEEHandle;

  iEvent.getByToken(recHitCollectionEBToken_, recHitCollectionEBHandle);
  iEvent.getByToken(recHitCollectionEEToken_, recHitCollectionEEHandle);

  std::unique_ptr<std::vector<T>> out = std::make_unique<std::vector<T>>();

  size_t nrObj = inHandle->size();
  std::array<std::vector<float>, EGEnergySysIndex::kNrSysErrs> results;
  for (auto& res : results)
    res.reserve(nrObj);

  const PhotonEnergyCalibrator::EventType evtType =
      iEvent.isRealData() ? PhotonEnergyCalibrator::EventType::DATA : PhotonEnergyCalibrator::EventType::MC;

  for (const auto& pho : *inHandle) {
    out->emplace_back(pho);

    if (semiDeterministicRng_)
      setSemiDetRandomSeed(iEvent, pho, nrObj, out->size());

    const EcalRecHitCollection* recHits =
        (pho.isEB()) ? recHitCollectionEBHandle.product() : recHitCollectionEEHandle.product();
    std::array<float, EGEnergySysIndex::kNrSysErrs> uncertainties =
        energyCorrector_.calibrate(out->back(), iEvent.id().run(), recHits, iEvent.streamID(), evtType);

    for (size_t index = 0; index < EGEnergySysIndex::kNrSysErrs; index++) {
      results[index].push_back(uncertainties[index]);
    }
  }

  auto fillAndStore = [&](auto handle) {
    for (const auto& mapToStore : valMapsToStore_) {
      fillAndStoreValueMap(iEvent, handle, results[mapToStore], EGEnergySysIndex::name(mapToStore));
    }
  };

  if (produceCalibratedObjs_) {
    fillAndStore(iEvent.put(std::move(out)));
  } else {
    fillAndStore(inHandle);
  }
}

//needs to be synced to CalibratedElectronProducers, want the same seed for a given SC
template <typename T>
void CalibratedPhotonProducerT<T>::setSemiDetRandomSeed(const edm::Event& iEvent,
                                                        const T& obj,
                                                        size_t nrObjs,
                                                        size_t objNr) {
  if (obj.superCluster().isNonnull()) {
    semiDeterministicRng_->SetSeed(egamma::getRandomSeedFromSC(iEvent, obj.superCluster()));
  } else {
    semiDeterministicRng_->SetSeed(egamma::getRandomSeedFromObj(iEvent, obj, nrObjs, objNr));
  }
}

using CalibratedPhotonProducer = CalibratedPhotonProducerT<reco::Photon>;
using CalibratedPatPhotonProducer = CalibratedPhotonProducerT<pat::Photon>;

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(CalibratedPhotonProducer);
DEFINE_FWK_MODULE(CalibratedPatPhotonProducer);