IsoValueMapProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PhysicsTools/NanoAOD
// Class:      IsoValueMapProducer
//
//
// Original Author:  Marco Peruzzi
//         Created:  Mon, 04 Sep 2017 22:43:53 GMT
//
//

// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"

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

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"

#include "RecoEgamma/EgammaTools/interface/EffectiveAreas.h"

#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "DataFormats/PatCandidates/interface/IsolatedTrack.h"

//
// class declaration
//

template <typename T>
class IsoValueMapProducer : public edm::global::EDProducer<> {
public:
  explicit IsoValueMapProducer(const edm::ParameterSet& iConfig)
      : src_(consumes<edm::View<T>>(iConfig.getParameter<edm::InputTag>("src"))),
        relative_(iConfig.getParameter<bool>("relative")) {
    if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
        typeid(T) == typeid(pat::IsolatedTrack)) {
      produces<edm::ValueMap<float>>("miniIsoChg");
      produces<edm::ValueMap<float>>("miniIsoAll");
      ea_miniiso_.reset(new EffectiveAreas((iConfig.getParameter<edm::FileInPath>("EAFile_MiniIso")).fullPath()));
      rho_miniiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_MiniIso"));
    }
    if ((typeid(T) == typeid(pat::Electron))) {
      produces<edm::ValueMap<float>>("PFIsoChg");
      produces<edm::ValueMap<float>>("PFIsoAll");
      produces<edm::ValueMap<float>>("PFIsoAll04");
      ea_pfiso_.reset(new EffectiveAreas((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso")).fullPath()));
      rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
    } else if ((typeid(T) == typeid(pat::Photon))) {
      produces<edm::ValueMap<float>>("PFIsoChg");
      produces<edm::ValueMap<float>>("PFIsoAll");
      ea_pfiso_chg_.reset(new EffectiveAreas((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Chg")).fullPath()));
      ea_pfiso_neu_.reset(new EffectiveAreas((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Neu")).fullPath()));
      ea_pfiso_pho_.reset(new EffectiveAreas((iConfig.getParameter<edm::FileInPath>("EAFile_PFIso_Pho")).fullPath()));
      rho_pfiso_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho_PFIso"));
    }
  }
  ~IsoValueMapProducer() override {}

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

  // ----------member data ---------------------------

  edm::EDGetTokenT<edm::View<T>> src_;
  bool relative_;
  edm::EDGetTokenT<double> rho_miniiso_;
  edm::EDGetTokenT<double> rho_pfiso_;
  std::unique_ptr<EffectiveAreas> ea_miniiso_;
  std::unique_ptr<EffectiveAreas> ea_pfiso_;
  std::unique_ptr<EffectiveAreas> ea_pfiso_chg_;
  std::unique_ptr<EffectiveAreas> ea_pfiso_neu_;
  std::unique_ptr<EffectiveAreas> ea_pfiso_pho_;
  float getEtaForEA(const T*) const;
  void doMiniIso(edm::Event&) const;
  void doPFIsoEle(edm::Event&) const;
  void doPFIsoPho(edm::Event&) const;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

template <typename T>
float IsoValueMapProducer<T>::getEtaForEA(const T* obj) const {
  return obj->eta();
}
template <>
float IsoValueMapProducer<pat::Electron>::getEtaForEA(const pat::Electron* el) const {
  return el->superCluster()->eta();
}
template <>
float IsoValueMapProducer<pat::Photon>::getEtaForEA(const pat::Photon* ph) const {
  return ph->superCluster()->eta();
}

template <typename T>
void IsoValueMapProducer<T>::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
  if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
      typeid(T) == typeid(pat::IsolatedTrack)) {
    doMiniIso(iEvent);
  };
  if ((typeid(T) == typeid(pat::Electron))) {
    doPFIsoEle(iEvent);
  }
  if ((typeid(T) == typeid(pat::Photon))) {
    doPFIsoPho(iEvent);
  }
}

template <typename T>
void IsoValueMapProducer<T>::doMiniIso(edm::Event& iEvent) const {
  edm::Handle<edm::View<T>> src;
  iEvent.getByToken(src_, src);
  edm::Handle<double> rho;
  iEvent.getByToken(rho_miniiso_, rho);

  unsigned int nInput = src->size();

  std::vector<float> miniIsoChg, miniIsoAll;
  miniIsoChg.reserve(nInput);
  miniIsoAll.reserve(nInput);

  for (const auto& obj : *src) {
    auto iso = obj.miniPFIsolation();
    auto chg = iso.chargedHadronIso();
    auto neu = iso.neutralHadronIso();
    auto pho = iso.photonIso();
    auto ea = ea_miniiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
    float R = 10.0 / std::min(std::max(obj.pt(), 50.0), 200.0);
    ea *= std::pow(R / 0.3, 2);
    float scale = relative_ ? 1.0 / obj.pt() : 1;
    miniIsoChg.push_back(scale * chg);
    miniIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - (*rho) * ea)));
  }

  std::unique_ptr<edm::ValueMap<float>> miniIsoChgV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerChg(*miniIsoChgV);
  fillerChg.insert(src, miniIsoChg.begin(), miniIsoChg.end());
  fillerChg.fill();
  std::unique_ptr<edm::ValueMap<float>> miniIsoAllV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerAll(*miniIsoAllV);
  fillerAll.insert(src, miniIsoAll.begin(), miniIsoAll.end());
  fillerAll.fill();

  iEvent.put(std::move(miniIsoChgV), "miniIsoChg");
  iEvent.put(std::move(miniIsoAllV), "miniIsoAll");
}

template <>
void IsoValueMapProducer<pat::Photon>::doMiniIso(edm::Event& iEvent) const {}

template <typename T>
void IsoValueMapProducer<T>::doPFIsoEle(edm::Event& iEvent) const {}

template <>
void IsoValueMapProducer<pat::Electron>::doPFIsoEle(edm::Event& iEvent) const {
  edm::Handle<edm::View<pat::Electron>> src;
  iEvent.getByToken(src_, src);
  edm::Handle<double> rho;
  iEvent.getByToken(rho_pfiso_, rho);

  unsigned int nInput = src->size();

  std::vector<float> PFIsoChg, PFIsoAll, PFIsoAll04;
  PFIsoChg.reserve(nInput);
  PFIsoAll.reserve(nInput);
  PFIsoAll04.reserve(nInput);

  for (const auto& obj : *src) {
    auto iso = obj.pfIsolationVariables();
    auto chg = iso.sumChargedHadronPt;
    auto neu = iso.sumNeutralHadronEt;
    auto pho = iso.sumPhotonEt;
    auto ea = ea_pfiso_->getEffectiveArea(fabs(getEtaForEA(&obj)));
    float scale = relative_ ? 1.0 / obj.pt() : 1;
    PFIsoChg.push_back(scale * chg);
    PFIsoAll.push_back(scale * (chg + std::max(0.0, neu + pho - (*rho) * ea)));
    PFIsoAll04.push_back(scale * (obj.chargedHadronIso() +
                                  std::max(0.0, obj.neutralHadronIso() + obj.photonIso() - (*rho) * ea * 16. / 9.)));
  }

  std::unique_ptr<edm::ValueMap<float>> PFIsoChgV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
  fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
  fillerChg.fill();
  std::unique_ptr<edm::ValueMap<float>> PFIsoAllV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
  fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
  fillerAll.fill();
  std::unique_ptr<edm::ValueMap<float>> PFIsoAll04V(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerAll04(*PFIsoAll04V);
  fillerAll04.insert(src, PFIsoAll04.begin(), PFIsoAll04.end());
  fillerAll04.fill();

  iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
  iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
  iEvent.put(std::move(PFIsoAll04V), "PFIsoAll04");
}

template <typename T>
void IsoValueMapProducer<T>::doPFIsoPho(edm::Event& iEvent) const {}

template <>
void IsoValueMapProducer<pat::Photon>::doPFIsoPho(edm::Event& iEvent) const {
  edm::Handle<edm::View<pat::Photon>> src;
  iEvent.getByToken(src_, src);
  edm::Handle<double> rho;
  iEvent.getByToken(rho_pfiso_, rho);

  unsigned int nInput = src->size();

  std::vector<float> PFIsoChg, PFIsoAll;
  PFIsoChg.reserve(nInput);
  PFIsoAll.reserve(nInput);

  for (const auto& obj : *src) {
    auto chg = obj.chargedHadronIso();
    auto neu = obj.neutralHadronIso();
    auto pho = obj.photonIso();
    auto ea_chg = ea_pfiso_chg_->getEffectiveArea(fabs(getEtaForEA(&obj)));
    auto ea_neu = ea_pfiso_neu_->getEffectiveArea(fabs(getEtaForEA(&obj)));
    auto ea_pho = ea_pfiso_pho_->getEffectiveArea(fabs(getEtaForEA(&obj)));
    float scale = relative_ ? 1.0 / obj.pt() : 1;
    PFIsoChg.push_back(scale * std::max(0.0, chg - (*rho) * ea_chg));
    PFIsoAll.push_back(PFIsoChg.back() +
                       scale * (std::max(0.0, neu - (*rho) * ea_neu) + std::max(0.0, pho - (*rho) * ea_pho)));
  }

  std::unique_ptr<edm::ValueMap<float>> PFIsoChgV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerChg(*PFIsoChgV);
  fillerChg.insert(src, PFIsoChg.begin(), PFIsoChg.end());
  fillerChg.fill();
  std::unique_ptr<edm::ValueMap<float>> PFIsoAllV(new edm::ValueMap<float>());
  edm::ValueMap<float>::Filler fillerAll(*PFIsoAllV);
  fillerAll.insert(src, PFIsoAll.begin(), PFIsoAll.end());
  fillerAll.fill();

  iEvent.put(std::move(PFIsoChgV), "PFIsoChg");
  iEvent.put(std::move(PFIsoAllV), "PFIsoAll");
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
template <typename T>
void IsoValueMapProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("src")->setComment("input physics object collection");
  desc.add<bool>("relative")->setComment("compute relative isolation instead of absolute one");
  if ((typeid(T) == typeid(pat::Muon)) || (typeid(T) == typeid(pat::Electron)) ||
      typeid(T) == typeid(pat::IsolatedTrack)) {
    desc.add<edm::FileInPath>("EAFile_MiniIso")
        ->setComment("txt file containing effective areas to be used for mini-isolation pileup subtraction");
    desc.add<edm::InputTag>("rho_MiniIso")
        ->setComment("rho to be used for effective-area based mini-isolation pileup subtraction");
  }
  if ((typeid(T) == typeid(pat::Electron))) {
    desc.add<edm::FileInPath>("EAFile_PFIso")
        ->setComment(
            "txt file containing effective areas to be used for PF-isolation pileup subtraction for electrons");
    desc.add<edm::InputTag>("rho_PFIso")
        ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for electrons");
  }
  if ((typeid(T) == typeid(pat::Photon))) {
    desc.add<edm::InputTag>("mapIsoChg")->setComment("input charged PF isolation calculated in VID for photons");
    desc.add<edm::InputTag>("mapIsoNeu")->setComment("input neutral PF isolation calculated in VID for photons");
    desc.add<edm::InputTag>("mapIsoPho")->setComment("input photon PF isolation calculated in VID for photons");
    desc.add<edm::FileInPath>("EAFile_PFIso_Chg")
        ->setComment(
            "txt file containing effective areas to be used for charged PF-isolation pileup subtraction for photons");
    desc.add<edm::FileInPath>("EAFile_PFIso_Neu")
        ->setComment(
            "txt file containing effective areas to be used for neutral PF-isolation pileup subtraction for photons");
    desc.add<edm::FileInPath>("EAFile_PFIso_Pho")
        ->setComment(
            "txt file containing effective areas to be used for photon PF-isolation pileup subtraction for photons");
    desc.add<edm::InputTag>("rho_PFIso")
        ->setComment("rho to be used for effective-area based PF-isolation pileup subtraction for photons");
  }
  std::string modname;
  if (typeid(T) == typeid(pat::Muon))
    modname += "Muon";
  else if (typeid(T) == typeid(pat::Electron))
    modname += "Ele";
  else if (typeid(T) == typeid(pat::Photon))
    modname += "Pho";
  else if (typeid(T) == typeid(pat::IsolatedTrack))
    modname += "IsoTrack";
  modname += "IsoValueMapProducer";
  descriptions.add(modname, desc);
}

typedef IsoValueMapProducer<pat::Muon> MuonIsoValueMapProducer;
typedef IsoValueMapProducer<pat::Electron> EleIsoValueMapProducer;
typedef IsoValueMapProducer<pat::Photon> PhoIsoValueMapProducer;
typedef IsoValueMapProducer<pat::IsolatedTrack> IsoTrackIsoValueMapProducer;

//define this as a plug-in
DEFINE_FWK_MODULE(MuonIsoValueMapProducer);
DEFINE_FWK_MODULE(EleIsoValueMapProducer);
DEFINE_FWK_MODULE(PhoIsoValueMapProducer);
DEFINE_FWK_MODULE(IsoTrackIsoValueMapProducer);