#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 &)

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

virtual	~ProducerBase ()

Private Types
typedef std::vector< T >	JetCollection

typedef edm::ValueMap< int >	JetSmearingFlags

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
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

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

Member Typedef Documentation

template<typename T , typename Textractor >

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

private

Definition at line 144 of file SmearedJetProducerT.h.

template<typename T , typename Textractor >

typedef edm::ValueMap<int> SmearedJetProducerT< T, Textractor >::JetSmearingFlags

private

Definition at line 145 of file SmearedJetProducerT.h.

Constructor & Destructor Documentation

template<typename T , typename Textractor >

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

inlineexplicit

Definition at line 149 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>();
     produces<JetSmearingFlags>("jetSmearingFlags");
 
   }

template<typename T , typename Textractor >

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

inline

Definition at line 202 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 211 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;
     }
 
     edm::Handle<JetCollection> jets;
     evt.getByLabel(src_, jets);
 
     int numJets = jets->size();
 
     std::auto_ptr<JetCollection> smearedJets(new JetCollection);
     std::vector<int> jetSmearingFlags_tmp(numJets);
 
     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();
       
       T rawJet(jet);
       rawJet.setP4(rawJetP4);
       double jetResolution = jetResolutionExtractor_(rawJet);
       double sigmaEn = jetResolution;
 
       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 ( TMath::Abs(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()));
       smearedJetEn = jet.energy() + (smearFactor - 1.)*dEn;
       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)
 
       double addSigmaEn = jetResolution*TMath::Sqrt(smearFactor*smearFactor - 1.);
       //smearedJetEn = jet.energy()*(1. + rnd_.Gaus(0., addSigmaEn)/TMath::Max(rawJetP4.E(), corrJetP4.E()));
       smearedJetEn = jet.energy() + rnd_.Gaus(0., addSigmaEn);
     }
       }
 
       // 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 << " multiplying 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);
 
       if ( isGenMatched ) jetSmearingFlags_tmp[jetIndex] = 1;
       else jetSmearingFlags_tmp[jetIndex] = 0;
     }
 
     std::auto_ptr<JetSmearingFlags> jetSmearingFlags(new JetSmearingFlags);
     JetSmearingFlags::Filler valueMapFiller(*jetSmearingFlags);
     valueMapFiller.insert(jets, jetSmearingFlags_tmp.begin(), jetSmearingFlags_tmp.end());
     valueMapFiller.fill();
 
 //--- add collection of "smeared" jets to the event
     evt.put(smearedJets);
     evt.put(jetSmearingFlags, "jetSmearingFlags");
   } 