MuonWithPFIsoProducer.cc

Go to the documentation of this file.

#include <memory>
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/EDProducer.h"

#include "DataFormats/Common/interface/View.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"

#include "MuonAnalysis/MomentumScaleCalibration/interface/MomentumScaleCorrector.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"

//
// class declaration
//
class MuonWithPFIsoProducer : public edm::EDProducer {
public:
  explicit MuonWithPFIsoProducer(const edm::ParameterSet&);
  ~MuonWithPFIsoProducer() override;

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

  edm::EDGetTokenT<edm::View<reco::Muon> > muonToken_;
  edm::EDGetTokenT<edm::View<reco::PFCandidate> > pfToken_;

  bool usePfMuonsOnly_;

  double trackIsoVeto_;
  double gammaIsoVeto_;
  double neutralHadronIsoVeto_;
};

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"

#include "DataFormats/GeometryVector/interface/VectorUtil.h"

MuonWithPFIsoProducer::MuonWithPFIsoProducer(const edm::ParameterSet& pset) {
  // What is being produced
  produces<std::vector<reco::Muon> >();

  // Muon collection
  muonToken_ =
      consumes<edm::View<reco::Muon> >(pset.getUntrackedParameter<edm::InputTag>("MuonTag", edm::InputTag("muons")));

  // PF candidate collection
  pfToken_ = consumes<edm::View<reco::PFCandidate> >(
      pset.getUntrackedParameter<edm::InputTag>("PFTag", edm::InputTag("particleFlow")));

  // Use only PF muons to get exact consistency with PfMET
  usePfMuonsOnly_ = pset.getUntrackedParameter<bool>("UsePfMuonsOnly", false);

  // Veto cone
  trackIsoVeto_ = pset.getUntrackedParameter<double>("TrackIsoVeto", 0.01);
  gammaIsoVeto_ = pset.getUntrackedParameter<double>("GammaIsoVeto", 0.07);
  neutralHadronIsoVeto_ = pset.getUntrackedParameter<double>("NeutralHadronIsoVeto", 0.1);
}

MuonWithPFIsoProducer::~MuonWithPFIsoProducer() {}

void MuonWithPFIsoProducer::beginJob() {}

void MuonWithPFIsoProducer::endJob() {}

void MuonWithPFIsoProducer::produce(edm::Event& ev, const edm::EventSetup& iSetup) {
  // Initialize pointer to new output muon collection
  std::unique_ptr<reco::MuonCollection> newmuons(new reco::MuonCollection);

  // Get Muon collection
  edm::Handle<edm::View<reco::Muon> > muonCollection;
  if (!ev.getByToken(muonToken_, muonCollection)) {
    edm::LogError("") << ">>> Muon collection does not exist !!!";
    ev.put(std::move(newmuons));
    return;
  }

  // Get PFCandidate collection
  edm::Handle<edm::View<reco::PFCandidate> > pfCollection;
  if (!ev.getByToken(pfToken_, pfCollection)) {
    edm::LogError("") << ">>> PFCandidate collection does not exist !!!";
    ev.put(std::move(newmuons));
    return;
  }

  // Loop over Pf candidates to find muons and collect deposits in veto,
  // dR<0.3 and dR<0.5 cones. Interpret "track" as charged particles (e,mu,
  // chraged hadrons). Interpret "em" as photons and also as electromagnetic
  // energy in HF. Interpret "had" as neutral hadrons and also as hadronic
  // energy in HF. Apply weights if requested at input level.
  // HO energies are not filled. Ditto for jet energies around the muon.
  unsigned int muonCollectionSize = muonCollection->size();
  unsigned int pfCollectionSize = pfCollection->size();
  for (unsigned int i = 0; i < muonCollectionSize; i++) {
    edm::RefToBase<reco::Muon> mu = muonCollection->refAt(i);

    // Ask for PfMuon consistency if requested
    bool muonFound = false;

    // Starting bycloning this muon
    reco::Muon* newmu = mu->clone();
    reco::TrackRef tk = mu->innerTrack();

    // Set isolations
    reco::MuonIsolation iso03;
    reco::MuonIsolation iso05;
    // Loop on all candidates
    for (unsigned int j = 0; j < pfCollectionSize; j++) {
      edm::RefToBase<reco::PFCandidate> pf = pfCollection->refAt(j);

      // Check the muon is in the PF collection when required
      bool thisIsTheMuon = false;
      if (tk.isNonnull() && pf->trackRef() == tk) {
        thisIsTheMuon = true;
        muonFound = true;
      }

      // Get dR. Nothing to add if dR>0.5
      double deltaR = Geom::deltaR(mu->momentum(), pf->momentum());
      if (deltaR > 0.5)
        continue;

      // Fill "tracker" components
      if (pf->particleId() == reco::PFCandidate::h || pf->particleId() == reco::PFCandidate::e ||
          pf->particleId() == reco::PFCandidate::mu) {
        if (deltaR < trackIsoVeto_ || thisIsTheMuon) {
          iso05.trackerVetoPt += pf->pt();
          iso03.trackerVetoPt += pf->pt();
        } else {
          iso05.sumPt += pf->pt();
          iso05.nTracks++;
          if (deltaR < 0.3) {
            iso03.sumPt += pf->pt();
            iso03.nTracks++;
          }
        }
        // Fill "em" components
      } else if (pf->particleId() == reco::PFCandidate::gamma || pf->particleId() == reco::PFCandidate::egamma_HF) {
        if (deltaR < gammaIsoVeto_) {
          iso05.emVetoEt += pf->pt();
          iso03.emVetoEt += pf->pt();
        } else {
          iso05.emEt += pf->pt();
          if (deltaR < 0.3)
            iso03.emEt += pf->pt();
        }
        // Fill "had" components
      } else if (pf->particleId() == reco::PFCandidate::h0 || pf->particleId() == reco::PFCandidate::h_HF) {
        if (deltaR < neutralHadronIsoVeto_) {
          iso05.hadVetoEt += pf->pt();
          iso03.hadVetoEt += pf->pt();
        } else {
          iso05.hadEt += pf->pt();
          if (deltaR < 0.3)
            iso03.hadEt += pf->pt();
        }
      }
    }

    // Do not take this muon (under explicit request) if it is not a PfMuon
    if (usePfMuonsOnly_ && (!muonFound))
      continue;

    // Set this isolation information in the new muon
    newmu->setIsolation(iso03, iso05);

    // Add new muon to output collection
    newmuons->push_back(*newmu);
  }

  // Add output collection to event
  ev.put(std::move(newmuons));
}

DEFINE_FWK_MODULE(MuonWithPFIsoProducer);