#include <SmearedJetProducerT.h>

Inheritance diagram for SmearedJetProducerT< T >:

Public Member Functions
virtual 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

Public Member Functions inherited from edm::stream::EDProducerBase
	EDProducerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDProducerBase ()

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector < edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, std::unordered_multimap< std::string, edm::ProductResolverIndex > const &iIndicies, std::string const &moduleLabel)

virtual	~ProducerBase () noexcept(false)

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector < ProductResolverIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

virtual	~EDConsumerBase () noexcept(false)

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

Static Public Member Functions inherited from edm::stream::EDProducerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (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

bool	m_use_txt_files

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

Public Types inherited from edm::stream::EDProducerBase
typedef EDProducerAdaptorBase	ModuleType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

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 99 of file SmearedJetProducerT.h.

Member Typedef Documentation

template<typename T >

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

private

Definition at line 101 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

template<typename T >

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

inlineexplicit

Definition at line 104 of file SmearedJetProducerT.h.

                                                                 :
             m_enabled(cfg.getParameter<bool>("enabled")),
             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_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 156 of file SmearedJetProducerT.h.

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

                                                                                  {
             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::uint32_t>("seed", 37428479);
             desc.add<bool>("skipGenMatching", false);
             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 >

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

inlineoverridevirtual

Implements edm::stream::EDProducerBase.

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;
 
             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);
                 }
             }
 
             const JetCollection& jets = *jets_collection;
 
             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::JetEta, jet.eta()}}, m_systematic_variation);
 
                 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));
         }