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

bool	hasAbilityToProduceInEndLumis () 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::ParameterSetDescription::add(), edm::ConfigurationDescriptions::addDefault(), edm::ParameterSetDescription::addNode(), edm::ParameterSetDescription::addUntracked(), 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));
   }