SmearedJetProducerT< T > Class Template Reference

#include <SmearedJetProducerT.h>

Inheritance diagram for SmearedJetProducerT< T >:

Public Member Functions
void	produce (edm::Event &event, const edm::EventSetup &setup) override
	SmearedJetProducerT (const edm::ParameterSet &cfg)
Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default
bool	hasAbilityToProduceInLumis () const final
bool	hasAbilityToProduceInRuns () const final

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &descriptions)

Private Types
using	JetCollection = std::vector< T >

Private Attributes
GreaterByPt< T >	jetPtComparator
bool	m_debug
bool	m_enabled
std::shared_ptr< pat::GenJetMatcher >	m_genJetMatcher
std::string	m_jets_algo
std::string	m_jets_algo_pt
edm::EDGetTokenT< JetCollection >	m_jets_token
int	m_nomVar
std::mt19937	m_random_generator
std::unique_ptr< JME::JetResolution >	m_resolution_from_file
edm::EDGetTokenT< double >	m_rho_token
std::unique_ptr< JME::JetResolutionScaleFactor >	m_scale_factor_from_file
Variation	m_systematic_variation
std::string	m_uncertaintySource
bool	m_use_txt_files
bool	m_useDeterministicSeed

Static Private Attributes
static const double	MIN_JET_ENERGY = 1e-2

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... >	CacheTypes
typedef CacheTypes::GlobalCache	GlobalCache
typedef AbilityChecker< T... >	HasAbility
typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache
typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext
typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache
typedef CacheTypes::RunCache	RunCache
typedef RunContextT< RunCache, GlobalCache >	RunContext
typedef CacheTypes::RunSummaryCache	RunSummaryCache

Detailed Description

template<typename T>
class SmearedJetProducerT< T >

Produce collection of "smeared" jets.

The aim of this correction is to account for the difference in jet energy resolution between Monte Carlo simulation and Data.

Author

Sébastien Brochet

Definition at line 97 of file SmearedJetProducerT.h.

Member Typedef Documentation

template<typename T >

using SmearedJetProducerT< T >::JetCollection = std::vector<T>

private

Definition at line 98 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

template<typename T >

SmearedJetProducerT< T >::SmearedJetProducerT ( const edm::ParameterSet &  cfg )

inlineexplicit

Definition at line 101 of file SmearedJetProducerT.h.

References looper::cfg, edm::errors::ConfigFileReadError, DOWN, Exception, edm::ParameterSet::exists(), edm::ParameterSet::getParameter(), NOMINAL, SurveyInfoScenario_cff::seed, patPFMETCorrections_cff::skipGenMatching, AlCaHLTBitMon_QueryRunRegistry::string, UP, and patPFMETCorrections_cff::variation.

      : m_enabled(cfg.getParameter<bool>("enabled")),
        m_useDeterministicSeed(cfg.getParameter<bool>("useDeterministicSeed")),
        m_debug(cfg.getUntrackedParameter<bool>("debug", false)) {
    m_jets_token = consumes<JetCollection>(cfg.getParameter<edm::InputTag>("src"));

    if (m_enabled) {
      m_rho_token = consumes<double>(cfg.getParameter<edm::InputTag>("rho"));

      m_use_txt_files = cfg.exists("resolutionFile") && cfg.exists("scaleFactorFile");

      if (m_use_txt_files) {
        std::string resolutionFile = cfg.getParameter<edm::FileInPath>("resolutionFile").fullPath();
        std::string scaleFactorFile = cfg.getParameter<edm::FileInPath>("scaleFactorFile").fullPath();

        m_resolution_from_file.reset(new JME::JetResolution(resolutionFile));
        m_scale_factor_from_file.reset(new JME::JetResolutionScaleFactor(scaleFactorFile));
      } else {
        m_jets_algo = cfg.getParameter<std::string>("algo");
        m_jets_algo_pt = cfg.getParameter<std::string>("algopt");
      }

      std::uint32_t seed = cfg.getParameter<std::uint32_t>("seed");
      m_random_generator = std::mt19937(seed);

      bool skipGenMatching = cfg.getParameter<bool>("skipGenMatching");
      if (!skipGenMatching)
        m_genJetMatcher = std::make_shared<pat::GenJetMatcher>(cfg, consumesCollector());

      std::int32_t variation = cfg.getParameter<std::int32_t>("variation");
      m_uncertaintySource = cfg.getParameter<std::string>("uncertaintySource");
      m_nomVar = 1;
      if (variation == 0)
        m_systematic_variation = Variation::NOMINAL;
      else if (variation == 1)
        m_systematic_variation = Variation::UP;
      else if (variation == -1)
        m_systematic_variation = Variation::DOWN;
      else if (variation == 101) {
        m_systematic_variation = Variation::NOMINAL;
        m_nomVar = 1;
      } else if (variation == -101) {
        m_systematic_variation = Variation::NOMINAL;
        m_nomVar = -1;
      } else
        throw edm::Exception(edm::errors::ConfigFileReadError,
                             "Invalid value for 'variation' parameter. Only -1, 0, 1 or 101, -101 are supported.");
    }

    produces<JetCollection>();
  }

Member Function Documentation

template<typename T >

static void SmearedJetProducerT< T >::fillDescriptions ( edm::ConfigurationDescriptions &  descriptions )

inlinestatic

Definition at line 153 of file SmearedJetProducerT.h.

References edm::ParameterSetDescription::add(), edm::ConfigurationDescriptions::addDefault(), edm::ParameterSetDescription::addNode(), edm::ParameterSetDescription::addUntracked(), pat::GenJetMatcher::fillDescriptions(), eostools::move(), source, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                           {
    edm::ParameterSetDescription desc;

    desc.add<edm::InputTag>("src");
    desc.add<bool>("enabled");
    desc.add<edm::InputTag>("rho");
    desc.add<std::int32_t>("variation", 0);
    desc.add<std::string>("uncertaintySource", "");
    desc.add<std::uint32_t>("seed", 37428479);
    desc.add<bool>("skipGenMatching", false);
    desc.add<bool>("useDeterministicSeed", true);
    desc.addUntracked<bool>("debug", false);

    auto source = (edm::ParameterDescription<std::string>("algo", true) and
                   edm::ParameterDescription<std::string>("algopt", true)) xor
                  (edm::ParameterDescription<edm::FileInPath>("resolutionFile", true) and
                   edm::ParameterDescription<edm::FileInPath>("scaleFactorFile", true));
    desc.addNode(std::move(source));

    pat::GenJetMatcher::fillDescriptions(desc);

    descriptions.addDefault(desc);
  }

template<typename T >

void SmearedJetProducerT< T >::produce ( edm::Event &  event, const edm::EventSetup &  setup  )

inlineoverride

Definition at line 177 of file SmearedJetProducerT.h.

References gather_cfg::cout, edmIntegrityCheck::d, reco::LeafCandidate::energy(), reco::LeafCandidate::eta(), JME::JetResolution::get(), JME::JetResolutionScaleFactor::get(), JME::JetResolution::getResolution(), JME::JetResolutionScaleFactor::getScaleFactor(), edm::EventBase::isRealData(), metsig::jet, JME::JetEta, JME::JetPt, fwrapper::jets, eostools::move(), reco::LeafCandidate::phi(), reco::LeafCandidate::pt(), fftjetproducer_cfi::resolution, JME::Rho, SurveyInfoScenario_cff::seed, jetUpdater_cfi::sort, and mathSSE::sqrt().

                                                                     {
    edm::Handle<JetCollection> jets_collection;
    event.getByToken(m_jets_token, jets_collection);

    // Disable the module when running on real data
    if (m_enabled && event.isRealData()) {
      m_enabled = false;
      m_genJetMatcher.reset();
    }

    edm::Handle<double> rho;
    if (m_enabled)
      event.getByToken(m_rho_token, rho);

    JME::JetResolution resolution;
    JME::JetResolutionScaleFactor resolution_sf;

    const JetCollection& jets = *jets_collection;

    if (m_enabled) {
      if (m_use_txt_files) {
        resolution = *m_resolution_from_file;
        resolution_sf = *m_scale_factor_from_file;
      } else {
        resolution = JME::JetResolution::get(setup, m_jets_algo_pt);
        resolution_sf = JME::JetResolutionScaleFactor::get(setup, m_jets_algo);
      }

      if (m_useDeterministicSeed) {
        unsigned int runNum_uint = static_cast<unsigned int>(event.id().run());
        unsigned int lumiNum_uint = static_cast<unsigned int>(event.id().luminosityBlock());
        unsigned int evNum_uint = static_cast<unsigned int>(event.id().event());
        unsigned int jet0eta = uint32_t(jets.empty() ? 0 : jets[0].eta() / 0.01);
        std::uint32_t seed = jet0eta + m_nomVar + (lumiNum_uint << 10) + (runNum_uint << 20) + evNum_uint;
        m_random_generator.seed(seed);
      }
    }

    if (m_genJetMatcher)
      m_genJetMatcher->getTokens(event);

    auto smearedJets = std::make_unique<JetCollection>();

    for (const auto& jet : jets) {
      if ((!m_enabled) || (jet.pt() == 0)) {
        // Module disabled or invalid p4. Simply copy the input jet.
        smearedJets->push_back(jet);

        continue;
      }

      double jet_resolution = resolution.getResolution(
          {{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}, {JME::Binning::Rho, *rho}});
      double jer_sf = resolution_sf.getScaleFactor({{JME::Binning::JetPt, jet.pt()}, {JME::Binning::JetEta, jet.eta()}},
                                                   m_systematic_variation,
                                                   m_uncertaintySource);
      if (m_debug) {
        std::cout << "jet:  pt: " << jet.pt() << "  eta: " << jet.eta() << "  phi: " << jet.phi()
                  << "  e: " << jet.energy() << std::endl;
        std::cout << "resolution: " << jet_resolution << std::endl;
        std::cout << "resolution scale factor: " << jer_sf << std::endl;
      }

      const reco::GenJet* genJet = nullptr;
      if (m_genJetMatcher)
        genJet = m_genJetMatcher->match(jet, jet.pt() * jet_resolution);

      double smearFactor = 1.;

      if (genJet) {
        /*
                     * Case 1: we have a "good" gen jet matched to the reco jet
                     */

        if (m_debug) {
          std::cout << "gen jet:  pt: " << genJet->pt() << "  eta: " << genJet->eta() << "  phi: " << genJet->phi()
                    << "  e: " << genJet->energy() << std::endl;
        }

        double dPt = jet.pt() - genJet->pt();
        smearFactor = 1 + m_nomVar * (jer_sf - 1.) * dPt / jet.pt();

      } else if (jer_sf > 1) {
        /*
                     * Case 2: we don't have a gen jet. Smear jet pt using a random gaussian variation
                     */

        double sigma = jet_resolution * std::sqrt(jer_sf * jer_sf - 1);
        if (m_debug) {
          std::cout << "gaussian width: " << sigma << std::endl;
        }

        std::normal_distribution<> d(0, sigma);
        smearFactor = 1. + m_nomVar * d(m_random_generator);
      } else if (m_debug) {
        std::cout << "Impossible to smear this jet" << std::endl;
      }

      if (jet.energy() * smearFactor < MIN_JET_ENERGY) {
        // Negative or too small smearFactor. We would change direction of the jet
        // and this is not what we want.
        // Recompute the smearing factor in order to have jet.energy() == MIN_JET_ENERGY
        double newSmearFactor = MIN_JET_ENERGY / jet.energy();
        if (m_debug) {
          std::cout << "The smearing factor (" << smearFactor << ") is either negative or too small. Fixing it to "
                    << newSmearFactor << " to avoid change of direction." << std::endl;
        }
        smearFactor = newSmearFactor;
      }

      T smearedJet = jet;
      smearedJet.scaleEnergy(smearFactor);

      if (m_debug) {
        std::cout << "smeared jet (" << smearFactor << "):  pt: " << smearedJet.pt() << "  eta: " << smearedJet.eta()
                  << "  phi: " << smearedJet.phi() << "  e: " << smearedJet.energy() << std::endl;
      }

      smearedJets->push_back(smearedJet);
    }

    // Sort jets by pt
    std::sort(smearedJets->begin(), smearedJets->end(), jetPtComparator);

    event.put(std::move(smearedJets));
  }

Member Data Documentation

template<typename T >

GreaterByPt<T> SmearedJetProducerT< T >::jetPtComparator

private

Definition at line 324 of file SmearedJetProducerT.h.

template<typename T >

bool SmearedJetProducerT< T >::m_debug

private

Definition at line 315 of file SmearedJetProducerT.h.

template<typename T >

bool SmearedJetProducerT< T >::m_enabled

private

Definition at line 309 of file SmearedJetProducerT.h.

template<typename T >

std::shared_ptr<pat::GenJetMatcher> SmearedJetProducerT< T >::m_genJetMatcher

private

Definition at line 316 of file SmearedJetProducerT.h.

template<typename T >

std::string SmearedJetProducerT< T >::m_jets_algo

private

Definition at line 311 of file SmearedJetProducerT.h.

template<typename T >

std::string SmearedJetProducerT< T >::m_jets_algo_pt

private

Definition at line 310 of file SmearedJetProducerT.h.

template<typename T >

edm::EDGetTokenT<JetCollection> SmearedJetProducerT< T >::m_jets_token

private

Definition at line 307 of file SmearedJetProducerT.h.

template<typename T >

int SmearedJetProducerT< T >::m_nomVar

private

Definition at line 326 of file SmearedJetProducerT.h.

template<typename T >

std::mt19937 SmearedJetProducerT< T >::m_random_generator

private

Definition at line 322 of file SmearedJetProducerT.h.

template<typename T >

std::unique_ptr<JME::JetResolution> SmearedJetProducerT< T >::m_resolution_from_file

private

Definition at line 319 of file SmearedJetProducerT.h.

template<typename T >

edm::EDGetTokenT<double> SmearedJetProducerT< T >::m_rho_token

private

Definition at line 308 of file SmearedJetProducerT.h.

template<typename T >

std::unique_ptr<JME::JetResolutionScaleFactor> SmearedJetProducerT< T >::m_scale_factor_from_file

private

Definition at line 320 of file SmearedJetProducerT.h.

template<typename T >

Variation SmearedJetProducerT< T >::m_systematic_variation

private

Definition at line 312 of file SmearedJetProducerT.h.

template<typename T >

std::string SmearedJetProducerT< T >::m_uncertaintySource

private

Definition at line 313 of file SmearedJetProducerT.h.

template<typename T >

bool SmearedJetProducerT< T >::m_use_txt_files

private

Definition at line 318 of file SmearedJetProducerT.h.

template<typename T >

bool SmearedJetProducerT< T >::m_useDeterministicSeed

private

Definition at line 314 of file SmearedJetProducerT.h.

template<typename T >

const double SmearedJetProducerT< T >::MIN_JET_ENERGY = 1e-2

staticprivate

Definition at line 305 of file SmearedJetProducerT.h.