DistortedMuonProducerFromDB.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 "FWCore/Framework/interface/Event.h"
 
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/MuonReco/interface/Muon.h"
 
#include "MuonAnalysis/MomentumScaleCalibration/interface/MomentumScaleCorrector.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/ResolutionFunction.h"
 
//
// class declaration
//
class DistortedMuonProducerFromDB : public edm::EDProducer {
public:
  explicit DistortedMuonProducerFromDB(const edm::ParameterSet&);
  ~DistortedMuonProducerFromDB() override;
 
private:
  void beginRun(const edm::Run&, const edm::EventSetup&) override;
  void produce(edm::Event&, const edm::EventSetup&) override;
  void endJob() override;
 
  edm::EDGetTokenT<edm::View<reco::Muon> > muonToken_;
 
  std::string dbScaleLabel_;
  std::string dbDataResolutionLabel_;
  std::string dbMCResolutionLabel_;
 
  std::unique_ptr<MomentumScaleCorrector> momCorrector_;
  std::unique_ptr<ResolutionFunction> momResolutionData_;
  std::unique_ptr<ResolutionFunction> momResolutionMC_;
};
 
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
 
#include <CLHEP/Random/RandGauss.h>
 
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "CondFormats/RecoMuonObjects/interface/MuScleFitDBobject.h"
#include "CondFormats/DataRecord/interface/MuScleFitDBobjectRcd.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/BaseFunction.h"
#include "MuonAnalysis/MomentumScaleCalibration/interface/Functions.h"
 
DistortedMuonProducerFromDB::DistortedMuonProducerFromDB(const edm::ParameterSet& pset) {
  // What is being produced
  produces<std::vector<reco::Muon> >();
 
  // Input products
  muonToken_ =
      consumes<edm::View<reco::Muon> >(pset.getUntrackedParameter<edm::InputTag>("MuonTag", edm::InputTag("muons")));
  dbScaleLabel_ = pset.getUntrackedParameter<std::string>("DBScaleLabel", "scale");
  dbDataResolutionLabel_ = pset.getUntrackedParameter<std::string>("DBDataResolutionLabel", "datareso");
  dbMCResolutionLabel_ = pset.getUntrackedParameter<std::string>("DBMCResolutionLabel", "mcreso");
}
 
DistortedMuonProducerFromDB::~DistortedMuonProducerFromDB() {}
 
void DistortedMuonProducerFromDB::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
  edm::ESHandle<MuScleFitDBobject> dbObject1;
  iSetup.get<MuScleFitDBobjectRcd>().get(dbScaleLabel_, dbObject1);
  momCorrector_ = std::make_unique<MomentumScaleCorrector>(dbObject1.product());
 
  LogTrace("") << ">>> Using database for momentum scale corrections !!";
 
  edm::ESHandle<MuScleFitDBobject> dbObject2;
  iSetup.get<MuScleFitDBobjectRcd>().get(dbDataResolutionLabel_, dbObject2);
  momResolutionData_ = std::make_unique<ResolutionFunction>(dbObject2.product());
 
  edm::ESHandle<MuScleFitDBobject> dbObject3;
  iSetup.get<MuScleFitDBobjectRcd>().get(dbMCResolutionLabel_, dbObject3);
  momResolutionMC_ = std::make_unique<ResolutionFunction>(dbObject3.product());
 
  LogTrace("") << ">>> Using database for momentum resolution corrections !!";
}
 
void DistortedMuonProducerFromDB::endJob() {}
 
void DistortedMuonProducerFromDB::produce(edm::Event& ev, const edm::EventSetup& iSetup) {
  if (ev.isRealData())
    return;
 
  // Muon collection
  edm::Handle<edm::View<reco::Muon> > muonCollection;
  if (!ev.getByToken(muonToken_, muonCollection)) {
    edm::LogError("") << ">>> Muon collection does not exist !!!";
    return;
  }
  unsigned int muonCollectionSize = muonCollection->size();
 
  std::unique_ptr<reco::MuonCollection> newmuons(new reco::MuonCollection);
 
  for (unsigned int i = 0; i < muonCollectionSize; i++) {
    edm::RefToBase<reco::Muon> mu = muonCollection->refAt(i);
 
    // Set shift
    double shift = (*momCorrector_)(*mu) - mu->pt();
    LogTrace("") << "\tmomentumScaleShift= " << shift << " [GeV]";
 
    // Set resolutions
    double sigma = pow(momResolutionData_->sigmaPt(*mu), 2) - pow(momResolutionMC_->sigmaPt(*mu), 2);
    if (sigma > 0.)
      sigma = sqrt(sigma);
    else
      sigma = 0.;
    LogTrace("") << "\tPt additional smearing= " << sigma << " [GeV]";
 
    // Gaussian Random number for smearing
    double rndg = CLHEP::RandGauss::shoot();
 
    // New muon
    double ptmu = mu->pt();
    ptmu += shift + sigma * rndg;
    reco::Muon* newmu = mu->clone();
    newmu->setP4(reco::Particle::PolarLorentzVector(ptmu, mu->eta(), mu->phi(), mu->mass()));
    newmuons->push_back(*newmu);
  }
 
  ev.put(std::move(newmuons));
}
 
DEFINE_FWK_MODULE(DistortedMuonProducerFromDB);