#include <SmearedJetProducerT.h>

Inheritance diagram for SmearedJetProducerT< T, Textractor >:

Public Member Functions
	SmearedJetProducerT (const edm::ParameterSet &cfg)

	~SmearedJetProducerT ()

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
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

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

virtual	~ProducerBase ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

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

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

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

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

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

Private Types
typedef std::vector< T >	JetCollection

Private Member Functions
virtual void	produce (edm::Event &evt, const edm::EventSetup &es)

Private Attributes
SmearedJetProducer_namespace::GenJetMatcherT < T >	genJetMatcher_

TFile *	inputFile_

double	jetCorrEtaMax_

Textractor	jetCorrExtractor_

edm::InputTag	jetCorrLabel_

edm::EDGetTokenT < reco::JetCorrector >	jetCorrToken_

SmearedJetProducer_namespace::JetResolutionExtractorT < T >	jetResolutionExtractor_

TH2 *	lut_

std::string	moduleLabel_

TRandom3	rnd_

double	shiftBy_

double	sigmaMaxGenJetMatch_

double	skipCorrJetPtThreshold_

StringCutObjectSelector< T > *	skipJetSelection_

double	skipRawJetPtThreshold_

double	smearBy_

edm::InputTag	src_

edm::EDGetTokenT< JetCollection >	srcToken_

int	verbosity_

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

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

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

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, typename Textractor>
class SmearedJetProducerT< T, Textractor >

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. The jet energy resolutions have been measured in QCD di-jet and gamma + jets events selected in 2010 data, as documented in the PAS JME-10-014.

Author: Christian Veelken, LLR

Definition at line 146 of file SmearedJetProducerT.h.

Member Typedef Documentation

template<typename T , typename Textractor >

typedef std::vector<T> SmearedJetProducerT< T, Textractor >::JetCollection

private

Definition at line 148 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

template<typename T , typename Textractor >

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

inlineexplicit

Definition at line 152 of file SmearedJetProducerT.h.

     : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
       genJetMatcher_(cfg, consumesCollector()),
       jetResolutionExtractor_(cfg.getParameter<edm::ParameterSet>("jetResolutions")),
       jetCorrLabel_(""),
       skipJetSelection_(0)
   {
     //std::cout << "<SmearedJetProducer::SmearedJetProducer>:" << std::endl;
     //std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
 
     src_ = cfg.getParameter<edm::InputTag>("src");
     srcToken_ = consumes<JetCollection>(src_);
 
     edm::FileInPath inputFileName = cfg.getParameter<edm::FileInPath>("inputFileName");
     std::string lutName = cfg.getParameter<std::string>("lutName");
     if (inputFileName.location() == edm::FileInPath::Unknown)
       throw cms::Exception("JetMETsmearInputProducer")
         << " Failed to find File = " << inputFileName << " !!\n";
 
     inputFile_ = new TFile(inputFileName.fullPath().data());
     lut_ = dynamic_cast<TH2*>(inputFile_->Get(lutName.data()));
     if ( !lut_ )
       throw cms::Exception("SmearedJetProducer")
         << " Failed to load LUT = " << lutName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";
 
     if ( cfg.exists("jetCorrLabel") ) {
       jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel");
       jetCorrToken_ = consumes<reco::JetCorrector>(jetCorrLabel_);
     }
 
     jetCorrEtaMax_ = ( cfg.exists("jetCorrEtaMax") ) ?
       cfg.getParameter<double>("jetCorrEtaMax") : 9.9;
 
     sigmaMaxGenJetMatch_ = cfg.getParameter<double>("sigmaMaxGenJetMatch");
 
     smearBy_ = ( cfg.exists("smearBy") ) ? cfg.getParameter<double>("smearBy") : 1.0;
     //std::cout << "smearBy = " << smearBy_ << std::endl;
 
     shiftBy_ = ( cfg.exists("shiftBy") ) ? cfg.getParameter<double>("shiftBy") : 0.;
     //std::cout << "shiftBy = " << shiftBy_ << std::endl;
 
     if ( cfg.exists("skipJetSelection") ) {
       std::string skipJetSelection_string = cfg.getParameter<std::string>("skipJetSelection");
       skipJetSelection_ = new StringCutObjectSelector<T>(skipJetSelection_string);
     }
 
     skipRawJetPtThreshold_  = ( cfg.exists("skipRawJetPtThreshold")  ) ?
       cfg.getParameter<double>("skipRawJetPtThreshold")  : 1.e-2;
     skipCorrJetPtThreshold_ = ( cfg.exists("skipCorrJetPtThreshold") ) ?
       cfg.getParameter<double>("skipCorrJetPtThreshold") : 1.e-2;
 
     verbosity_ = ( cfg.exists("verbosity") ) ?
       cfg.getParameter<int>("verbosity") : 0;
 
     produces<JetCollection>();
   }

template<typename T , typename Textractor >

SmearedJetProducerT< T, Textractor >::~SmearedJetProducerT ( )

inline

Definition at line 208 of file SmearedJetProducerT.h.

References SmearedJetProducerT< T, Textractor >::inputFile_, SmearedJetProducerT< T, Textractor >::lut_, and SmearedJetProducerT< T, Textractor >::skipJetSelection_.

   {
     delete skipJetSelection_;
     delete inputFile_;
     delete lut_;
   }

Member Function Documentation

template<typename T , typename Textractor >

virtual void SmearedJetProducerT< T, Textractor >::produce	(	edm::Event &	evt,
		const edm::EventSetup &	es
	)

inlineprivatevirtual

Implements edm::stream::EDProducerBase.

Definition at line 217 of file SmearedJetProducerT.h.

   {
     if ( verbosity_ ) {
       std::cout << "<SmearedJetProducerT::produce>:" << std::endl;
       std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
       std::cout << " src = " << src_.label() << std::endl;
     }
 
     std::auto_ptr<JetCollection> smearedJets(new JetCollection);
 
     edm::Handle<JetCollection> jets;
     evt.getByToken(srcToken_, jets);
 
     int numJets = jets->size();
     for ( int jetIndex = 0; jetIndex < numJets; ++jetIndex ) {
       const T& jet = jets->at(jetIndex);
 
       static const SmearedJetProducer_namespace::RawJetExtractorT<T> rawJetExtractor;
       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
       if ( verbosity_ ) {
     std::cout << "rawJet: Pt = " << rawJetP4.pt() << ", eta = " << rawJetP4.eta() << ", phi = " << rawJetP4.phi() << std::endl;
       }
 
       reco::Candidate::LorentzVector corrJetP4 = jet.p4();
       if ( !jetCorrLabel_.label().empty() ) {
         edm::Handle<reco::JetCorrector> jetCorr;
         evt.getByToken(jetCorrToken_, jetCorr);
     corrJetP4 =
     std::is_base_of<PATJetCorrExtractor, Textractor>::value ?
       jetCorrExtractor_(jet, jetCorrLabel_.label(), jetCorrEtaMax_, &rawJetP4) :
       jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
       }
       if ( verbosity_ ) {
     std::cout << "corrJet: Pt = " << corrJetP4.pt() << ", eta = " << corrJetP4.eta() << ", phi = " << corrJetP4.phi() << std::endl;
       }
 
       double smearFactor = 1.;
       double x = std::abs(corrJetP4.eta());
       double y = corrJetP4.pt();
       if ( x > lut_->GetXaxis()->GetXmin() && x < lut_->GetXaxis()->GetXmax() &&
        y > lut_->GetYaxis()->GetXmin() && y < lut_->GetYaxis()->GetXmax() ) {
     int binIndex = lut_->FindBin(x, y);
 
     if ( smearBy_ > 0. ) smearFactor += smearBy_*(lut_->GetBinContent(binIndex) - 1.);
     double smearFactorErr = lut_->GetBinError(binIndex);
     if ( verbosity_ ) std::cout << "smearFactor = " << smearFactor << " +/- " << smearFactorErr << std::endl;
 
     if ( shiftBy_ != 0. ) {
       smearFactor += (shiftBy_*smearFactorErr);
       if ( verbosity_ ) std::cout << "smearFactor(shifted) = " << smearFactor << std::endl;
     }
       }
 
       double smearedJetEn = jet.energy();
       double sigmaEn = jetResolutionExtractor_(jet)*sqrt(smearFactor*smearFactor - 1.);
       const reco::GenJet* genJet = genJetMatcher_(jet, &evt);
       bool isGenMatched = false;
       if ( genJet ) {
     if ( verbosity_ ) {
       std::cout << "genJet: Pt = " << genJet->pt() << ", eta = " << genJet->eta() << ", phi = " << genJet->phi() << std::endl;
     }
     double dEn = corrJetP4.E() - genJet->energy();
     if ( dEn < (sigmaMaxGenJetMatch_*sigmaEn) ) {
 //--- case 1: reconstructed jet matched to generator level jet,
 //            smear difference between reconstructed and "true" jet energy
 
       if ( verbosity_ ) {
         std::cout << " successfully matched to genJet" << std::endl;
         std::cout << "corrJetEn = " << corrJetP4.E() << ", genJetEn = " << genJet->energy() << " --> dEn = " << dEn << std::endl;
       }
 
       smearedJetEn = jet.energy()*(1. + (smearFactor - 1.)*dEn/std::max( rawJetP4.E(), corrJetP4.E()));
       isGenMatched = true;
     }
       }
       if ( !isGenMatched ) {
 //--- case 2: reconstructed jet **not** matched to generator level jet,
 //            smear jet energy using MC resolution functions implemented in PFMEt significance algorithm (CMS AN-10/400)
 
     if ( verbosity_ ) {
       std::cout << " not matched to genJet" << std::endl;
       std::cout << "corrJetEn = " << corrJetP4.E() << ", sigmaEn = " << sigmaEn << std::endl;
     }
 
     if ( smearFactor > 1. ) {
       // CV: MC resolution already accounted for in reconstructed jet,
       //     add additional Gaussian smearing of width = sqrt(smearFactor^2 - 1)
       //     to account for Data/MC **difference** in jet resolutions.
       //     Take maximum(rawJetEn, corrJetEn) to avoid pathological cases
       //    (e.g. corrJetEn << rawJetEn, due to L1Fastjet corrections)
 
       smearedJetEn = jet.energy()*(1. + rnd_.Gaus(0., sigmaEn)/std::max( rawJetP4.E(), corrJetP4.E()) );
     }
       }
 
       // CV: keep minimum jet energy, in order not to loose direction information
       const double minJetEn = 1.e-2;
       if ( smearedJetEn < minJetEn ) smearedJetEn = minJetEn;
 
       // CV: skip smearing in case either "raw" or "corrected" jet energy is very low
       //     or jet passes selection configurable via python
       //    (allows for protection against "pathological cases",
       //     cf. PhysicsTools/PatUtils/python/tools/metUncertaintyTools.py)
       reco::Candidate::LorentzVector smearedJetP4 = jet.p4();
       if ( !((skipJetSelection_ && (*skipJetSelection_)(jet)) ||
          rawJetP4.pt()  < skipRawJetPtThreshold_          ||
          corrJetP4.pt() < skipCorrJetPtThreshold_         ) ) {
     if ( verbosity_ ) {
       std::cout << " smearing jetP4 by factor = " << (smearedJetEn/jet.energy()) << " --> smearedJetEn = " << smearedJetEn << std::endl;
     }
     smearedJetP4 *= (smearedJetEn/jet.energy());
       }
 
       if ( verbosity_ ) {
     std::cout << "smearedJet: Pt = " << smearedJetP4.pt() << ", eta = " << smearedJetP4.eta() << ", phi = " << smearedJetP4.phi() << std::endl;
     std::cout << " dPt = " << (smearedJetP4.pt() - jet.pt())
           << " (Px = " << (smearedJetP4.px() - jet.px()) << ", Py = " << (smearedJetP4.py() - jet.py()) << ")" << std::endl;
       }
 
       T smearedJet = (jet);
       smearedJet.setP4(smearedJetP4);
 
       smearedJets->push_back(smearedJet);
     }
 
 //--- add collection of "smeared" jets to the event
     evt.put(smearedJets);
   }