#include <SmearedJetProducerT.h>

Inheritance diagram for SmearedJetProducerT< T, Textractor >:

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

	~SmearedJetProducerT ()

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

virtual	~EDProducer ()

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

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

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

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

void	labelsForToken (EDGetToken iToken, Labels &oLabels) 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_

std::string	jetCorrLabel_

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_

int	verbosity_

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

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

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

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

Member Typedef Documentation

template<typename T , typename Textractor >

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

private

Definition at line 143 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

template<typename T , typename Textractor >

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

inlineexplicit

Definition at line 147 of file SmearedJetProducerT.h.

     : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
       genJetMatcher_(cfg),
       jetResolutionExtractor_(cfg.getParameter<edm::ParameterSet>("jetResolutions")),
       skipJetSelection_(0)
   {
     //std::cout << "<SmearedJetProducer::SmearedJetProducer>:" << std::endl;
     //std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
 
     src_ = cfg.getParameter<edm::InputTag>("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";
 
     jetCorrLabel_ = ( cfg.exists("jetCorrLabel") ) ?
       cfg.getParameter<std::string>("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 198 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::EDProducer.

Definition at line 207 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.getByLabel(src_, jets);
 
     int numJets = jets->size();
     for ( int jetIndex = 0; jetIndex < numJets; ++jetIndex ) {
       const T& jet = jets->at(jetIndex);
       
       static 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_ != "" ) corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_, &evt, &es, jetCorrEtaMax_, &rawJetP4);
       if ( verbosity_ ) {
     std::cout << "corrJet: Pt = " << corrJetP4.pt() << ", eta = " << corrJetP4.eta() << ", phi = " << corrJetP4.phi() << std::endl;
       }
 
       double smearFactor = 1.;      
       double x = TMath::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)*TMath::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/TMath::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)/TMath::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);
   } 