#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	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () 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 constexpr 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

◆ JetCollection

template<typename T >

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

private

Definition at line 98 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

◆ SmearedJetProducerT()

template<typename T >

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

inlineexplicit

Definition at line 101 of file SmearedJetProducerT.h.

       : 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

◆ fillDescriptions()

template<typename T >

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

inlinestatic

Definition at line 153 of file SmearedJetProducerT.h.

                                                                            {
     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);
   }

References edm::ConfigurationDescriptions::addDefault(), submitPVResolutionJobs::desc, pat::GenJetMatcher::fillDescriptions(), eostools::move(), source, and AlCaHLTBitMon_QueryRunRegistry::string.

◆ produce()

template<typename T >

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

inlineoverride

Definition at line 177 of file SmearedJetProducerT.h.

                                                                      {
     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));
   }