#include <SmearedJetProducerT.h>

Inheritance diagram for SmearedJetProducerT< T, Textractor >:

Public Member Functions
	SmearedJetProducerT (const edm::ParameterSet &cfg)
	~SmearedJetProducerT ()
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_

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

Version:

Revision:: 1.9

Id:: SmearedJetProducerT.h,v 1.9 2012/08/31 09:58:44 veelken Exp

Definition at line 143 of file SmearedJetProducerT.h.

Member Typedef Documentation

template<typename T , typename Textractor >

typedef std::vector<T> SmearedJetProducerT< T, Textractor >::JetCollection [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 ) [inline, explicit]

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.isLocal() ) 
      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 200 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
	)		`[inline, private, virtual]`

Implements edm::EDProducer.

Definition at line 209 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);
  }