CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/PhysicsTools/PatUtils/interface/SmearedJetProducerT.h

Go to the documentation of this file.
00001 #ifndef PhysicsTools_PatUtils_SmearedJetProducerT_h
00002 #define PhysicsTools_PatUtils_SmearedJetProducerT_h
00003 
00018 #include "FWCore/Framework/interface/Frameworkfwd.h"
00019 #include "FWCore/Framework/interface/EDProducer.h"
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00022 #include "FWCore/Utilities/interface/InputTag.h"
00023 
00024 #include "FWCore/ParameterSet/interface/FileInPath.h"
00025 #include "FWCore/Utilities/interface/Exception.h"
00026 
00027 #include "DataFormats/PatCandidates/interface/Jet.h"
00028 #include "DataFormats/METReco/interface/CorrMETData.h"
00029 #include "DataFormats/Common/interface/Handle.h"
00030 
00031 #include "DataFormats/JetReco/interface/GenJet.h"
00032 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00033 #include "DataFormats/Math/interface/deltaR.h"
00034 
00035 #include "JetMETCorrections/Type1MET/interface/JetCorrExtractorT.h"
00036 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
00037 
00038 #include <TFile.h>
00039 #include <TFormula.h>
00040 #include <TH2.h>
00041 #include <TMath.h>
00042 #include <TRandom3.h>
00043 #include <TString.h>
00044 
00045 #include <iostream>
00046 #include <iomanip>
00047 
00048 namespace SmearedJetProducer_namespace
00049 {
00050   template <typename T>
00051   class GenJetMatcherT
00052   {
00053     public:
00054 
00055      GenJetMatcherT(const edm::ParameterSet& cfg) 
00056        : srcGenJets_(cfg.getParameter<edm::InputTag>("srcGenJets")),
00057          dRmaxGenJetMatch_(0)
00058      {
00059        TString dRmaxGenJetMatch_formula = cfg.getParameter<std::string>("dRmaxGenJetMatch").data();
00060        dRmaxGenJetMatch_formula.ReplaceAll("genJetPt", "x");
00061        dRmaxGenJetMatch_ = new TFormula("dRmaxGenJetMatch", dRmaxGenJetMatch_formula.Data());
00062      }
00063      ~GenJetMatcherT() 
00064      {
00065        delete dRmaxGenJetMatch_;
00066      }
00067 
00068      const reco::GenJet* operator()(const T& jet, edm::Event* evt = 0) const
00069      {
00070        assert(evt);
00071        
00072        edm::Handle<reco::GenJetCollection> genJets;
00073        evt->getByLabel(srcGenJets_, genJets);
00074 
00075        const reco::GenJet* retVal = 0;
00076 
00077        double dRbestMatch = 1.e+6;
00078        for ( reco::GenJetCollection::const_iterator genJet = genJets->begin();
00079              genJet != genJets->end(); ++genJet ) {
00080          double dRmax = dRmaxGenJetMatch_->Eval(genJet->pt());
00081          //std::cout << "genJetPt = " << genJet->pt() << ": dRmax = " << dRmax << std::endl;
00082          double dR = deltaR(jet.p4(), genJet->p4());     
00083          if ( dR < dRbestMatch && dR < dRmax ) {
00084            retVal = &(*genJet);
00085            dRbestMatch = dR;
00086          }
00087        }
00088 
00089        return retVal;
00090      }
00091 
00092     private:
00093 
00094 //--- configuration parameter
00095      edm::InputTag srcGenJets_;
00096 
00097      TFormula* dRmaxGenJetMatch_;
00098   };
00099 
00100   template <typename T>
00101   class JetResolutionExtractorT
00102   {
00103     public:
00104 
00105      JetResolutionExtractorT(const edm::ParameterSet&) {}
00106      ~JetResolutionExtractorT() {}
00107 
00108      double operator()(const T&) const
00109      {
00110        throw cms::Exception("SmearedJetProducer::produce")
00111          << " Jets of type other than PF not supported yet !!\n";       
00112      }
00113   };
00114 
00115   template <typename T>
00116   class RawJetExtractorT // this template is neccessary to support pat::Jets
00117                          // (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
00118   {
00119     public:
00120 
00121      reco::Candidate::LorentzVector operator()(const T& jet) const 
00122      {        
00123        return jet.p4();
00124      } 
00125   };
00126 
00127   template <>
00128   class RawJetExtractorT<pat::Jet>
00129   {
00130     public:
00131 
00132      reco::Candidate::LorentzVector operator()(const pat::Jet& jet) const 
00133      { 
00134        if ( jet.jecSetsAvailable() ) return jet.correctedP4("Uncorrected");
00135        else return jet.p4();
00136      } 
00137   };
00138 }
00139 
00140 template <typename T, typename Textractor>
00141 class SmearedJetProducerT : public edm::EDProducer 
00142 {
00143   typedef std::vector<T> JetCollection;
00144 
00145  public:
00146 
00147   explicit SmearedJetProducerT(const edm::ParameterSet& cfg)
00148     : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
00149       genJetMatcher_(cfg),
00150       jetResolutionExtractor_(cfg.getParameter<edm::ParameterSet>("jetResolutions")),
00151       skipJetSelection_(0)
00152   {
00153     //std::cout << "<SmearedJetProducer::SmearedJetProducer>:" << std::endl;
00154     //std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
00155 
00156     src_ = cfg.getParameter<edm::InputTag>("src");
00157 
00158     edm::FileInPath inputFileName = cfg.getParameter<edm::FileInPath>("inputFileName");
00159     std::string lutName = cfg.getParameter<std::string>("lutName");
00160     if (inputFileName.location() == edm::FileInPath::Unknown)
00161       throw cms::Exception("JetMETsmearInputProducer") 
00162         << " Failed to find File = " << inputFileName << " !!\n";
00163 
00164     inputFile_ = new TFile(inputFileName.fullPath().data());
00165     lut_ = dynamic_cast<TH2*>(inputFile_->Get(lutName.data()));
00166     if ( !lut_ ) 
00167       throw cms::Exception("SmearedJetProducer") 
00168         << " Failed to load LUT = " << lutName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";
00169 
00170     jetCorrLabel_ = ( cfg.exists("jetCorrLabel") ) ?
00171       cfg.getParameter<std::string>("jetCorrLabel") : "";
00172     jetCorrEtaMax_ = ( cfg.exists("jetCorrEtaMax") ) ?
00173       cfg.getParameter<double>("jetCorrEtaMax") : 9.9;
00174 
00175     sigmaMaxGenJetMatch_ = cfg.getParameter<double>("sigmaMaxGenJetMatch");
00176 
00177     smearBy_ = ( cfg.exists("smearBy") ) ? cfg.getParameter<double>("smearBy") : 1.0;
00178     //std::cout << "smearBy = " << smearBy_ << std::endl;
00179 
00180     shiftBy_ = ( cfg.exists("shiftBy") ) ? cfg.getParameter<double>("shiftBy") : 0.;
00181     //std::cout << "shiftBy = " << shiftBy_ << std::endl;
00182 
00183     if ( cfg.exists("skipJetSelection") ) {
00184       std::string skipJetSelection_string = cfg.getParameter<std::string>("skipJetSelection");
00185       skipJetSelection_ = new StringCutObjectSelector<T>(skipJetSelection_string);
00186     }
00187 
00188     skipRawJetPtThreshold_  = ( cfg.exists("skipRawJetPtThreshold")  ) ? 
00189       cfg.getParameter<double>("skipRawJetPtThreshold")  : 1.e-2;
00190     skipCorrJetPtThreshold_ = ( cfg.exists("skipCorrJetPtThreshold") ) ? 
00191       cfg.getParameter<double>("skipCorrJetPtThreshold") : 1.e-2;
00192 
00193     verbosity_ = ( cfg.exists("verbosity") ) ?
00194       cfg.getParameter<int>("verbosity") : 0;
00195 
00196     produces<JetCollection>();
00197   }
00198   ~SmearedJetProducerT()
00199   {
00200     delete skipJetSelection_;
00201     delete inputFile_;
00202     delete lut_;
00203   }
00204     
00205  private:
00206 
00207   virtual void produce(edm::Event& evt, const edm::EventSetup& es)
00208   {
00209     if ( verbosity_ ) {
00210       std::cout << "<SmearedJetProducerT::produce>:" << std::endl;
00211       std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
00212       std::cout << " src = " << src_.label() << std::endl;
00213     }
00214 
00215     std::auto_ptr<JetCollection> smearedJets(new JetCollection);
00216     
00217     edm::Handle<JetCollection> jets;
00218     evt.getByLabel(src_, jets);
00219 
00220     int numJets = jets->size();
00221     for ( int jetIndex = 0; jetIndex < numJets; ++jetIndex ) {
00222       const T& jet = jets->at(jetIndex);
00223       
00224       static SmearedJetProducer_namespace::RawJetExtractorT<T> rawJetExtractor;
00225       reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
00226       if ( verbosity_ ) {
00227         std::cout << "rawJet: Pt = " << rawJetP4.pt() << ", eta = " << rawJetP4.eta() << ", phi = " << rawJetP4.phi() << std::endl;
00228       }
00229 
00230       reco::Candidate::LorentzVector corrJetP4 = jet.p4();
00231       if ( jetCorrLabel_ != "" ) corrJetP4 = jetCorrExtractor_(jet, jetCorrLabel_, &evt, &es, jetCorrEtaMax_, &rawJetP4);
00232       if ( verbosity_ ) {
00233         std::cout << "corrJet: Pt = " << corrJetP4.pt() << ", eta = " << corrJetP4.eta() << ", phi = " << corrJetP4.phi() << std::endl;
00234       }
00235 
00236       double smearFactor = 1.;      
00237       double x = TMath::Abs(corrJetP4.eta());
00238       double y = corrJetP4.pt();
00239       if ( x > lut_->GetXaxis()->GetXmin() && x < lut_->GetXaxis()->GetXmax() &&
00240            y > lut_->GetYaxis()->GetXmin() && y < lut_->GetYaxis()->GetXmax() ) {
00241         int binIndex = lut_->FindBin(x, y);
00242         
00243         if ( smearBy_ > 0. ) smearFactor += smearBy_*(lut_->GetBinContent(binIndex) - 1.);
00244         double smearFactorErr = lut_->GetBinError(binIndex);
00245         if ( verbosity_ ) std::cout << "smearFactor = " << smearFactor << " +/- " << smearFactorErr << std::endl;
00246 
00247         if ( shiftBy_ != 0. ) {
00248           smearFactor += (shiftBy_*smearFactorErr);
00249           if ( verbosity_ ) std::cout << "smearFactor(shifted) = " << smearFactor << std::endl;
00250         }
00251       }
00252 
00253       double smearedJetEn = jet.energy();
00254       double sigmaEn = jetResolutionExtractor_(jet)*TMath::Sqrt(smearFactor*smearFactor - 1.);
00255       const reco::GenJet* genJet = genJetMatcher_(jet, &evt);
00256       bool isGenMatched = false;
00257       if ( genJet ) {
00258         if ( verbosity_ ) {
00259           std::cout << "genJet: Pt = " << genJet->pt() << ", eta = " << genJet->eta() << ", phi = " << genJet->phi() << std::endl;
00260         }
00261         double dEn = corrJetP4.E() - genJet->energy();
00262         if ( dEn < (sigmaMaxGenJetMatch_*sigmaEn) ) {
00263 //--- case 1: reconstructed jet matched to generator level jet, 
00264 //            smear difference between reconstructed and "true" jet energy
00265 
00266           if ( verbosity_ ) {
00267             std::cout << " successfully matched to genJet" << std::endl;        
00268             std::cout << "corrJetEn = " << corrJetP4.E() << ", genJetEn = " << genJet->energy() << " --> dEn = " << dEn << std::endl;
00269           }
00270 
00271           smearedJetEn = jet.energy()*(1. + (smearFactor - 1.)*dEn/TMath::Max(rawJetP4.E(), corrJetP4.E()));
00272           isGenMatched = true;
00273         }
00274       }
00275       if ( !isGenMatched ) {
00276 //--- case 2: reconstructed jet **not** matched to generator level jet, 
00277 //            smear jet energy using MC resolution functions implemented in PFMEt significance algorithm (CMS AN-10/400)
00278 
00279         if ( verbosity_ ) {
00280           std::cout << " not matched to genJet" << std::endl;
00281           std::cout << "corrJetEn = " << corrJetP4.E() << ", sigmaEn = " << sigmaEn << std::endl;
00282         }
00283 
00284         if ( smearFactor > 1. ) {
00285           // CV: MC resolution already accounted for in reconstructed jet,
00286           //     add additional Gaussian smearing of width = sqrt(smearFactor^2 - 1) 
00287           //     to account for Data/MC **difference** in jet resolutions.
00288           //     Take maximum(rawJetEn, corrJetEn) to avoid pathological cases
00289           //    (e.g. corrJetEn << rawJetEn, due to L1Fastjet corrections)
00290 
00291           smearedJetEn = jet.energy()*(1. + rnd_.Gaus(0., sigmaEn)/TMath::Max(rawJetP4.E(), corrJetP4.E()));
00292         }
00293       }
00294 
00295       // CV: keep minimum jet energy, in order not to loose direction information
00296       const double minJetEn = 1.e-2;
00297       if ( smearedJetEn < minJetEn ) smearedJetEn = minJetEn;
00298 
00299       // CV: skip smearing in case either "raw" or "corrected" jet energy is very low
00300       //     or jet passes selection configurable via python
00301       //    (allows for protection against "pathological cases",
00302       //     cf. PhysicsTools/PatUtils/python/tools/metUncertaintyTools.py)
00303       reco::Candidate::LorentzVector smearedJetP4 = jet.p4();
00304       if ( !((skipJetSelection_ && (*skipJetSelection_)(jet)) ||
00305              rawJetP4.pt()  < skipRawJetPtThreshold_          ||
00306              corrJetP4.pt() < skipCorrJetPtThreshold_         ) ) {
00307         if ( verbosity_ ) {
00308           std::cout << " smearing jetP4 by factor = " << (smearedJetEn/jet.energy()) << " --> smearedJetEn = " << smearedJetEn << std::endl;
00309         }
00310         smearedJetP4 *= (smearedJetEn/jet.energy());
00311       }
00312           
00313       if ( verbosity_ ) {
00314         std::cout << "smearedJet: Pt = " << smearedJetP4.pt() << ", eta = " << smearedJetP4.eta() << ", phi = " << smearedJetP4.phi() << std::endl;
00315         std::cout << " dPt = " << (smearedJetP4.pt() - jet.pt()) 
00316                   << " (Px = " << (smearedJetP4.px() - jet.px()) << ", Py = " << (smearedJetP4.py() - jet.py()) << ")" << std::endl;
00317       }
00318       
00319       T smearedJet = (jet);
00320       smearedJet.setP4(smearedJetP4);
00321       
00322       smearedJets->push_back(smearedJet);
00323     }
00324 
00325 //--- add collection of "smeared" jets to the event
00326     evt.put(smearedJets);
00327   } 
00328 
00329   std::string moduleLabel_;
00330       
00331   SmearedJetProducer_namespace::GenJetMatcherT<T> genJetMatcher_;
00332 
00333 //--- configuration parameters
00334  
00335   // collection of pat::Jets (with L2L3/L2L3Residual corrections applied)
00336   edm::InputTag src_;
00337 
00338   TFile* inputFile_;
00339   TH2* lut_;
00340 
00341   SmearedJetProducer_namespace::JetResolutionExtractorT<T> jetResolutionExtractor_;
00342   TRandom3 rnd_;
00343 
00344   std::string jetCorrLabel_; // e.g. 'ak5PFJetL1FastL2L3' (reco::PFJets) / '' (pat::Jets)
00345   double jetCorrEtaMax_; // do not use JEC factors for |eta| above this threshold (recommended default = 4.7),
00346                          // in order to work around problem with CMSSW_4_2_x JEC factors at high eta,
00347                          // reported in
00348                          //  https://hypernews.cern.ch/HyperNews/CMS/get/jes/270.html
00349                          //  https://hypernews.cern.ch/HyperNews/CMS/get/JetMET/1259/1.html
00350   Textractor jetCorrExtractor_;
00351 
00352   double sigmaMaxGenJetMatch_; // maximum difference between energy of reconstructed jet and matched generator level jet
00353                                // (if the difference between reconstructed and generated jet energy exceeds this threshold,
00354                                //  the jet is considered to have substantial pile-up contributions are is considered to be unmatched)
00355 
00356   double smearBy_; // option to "smear" jet energy by N standard-deviations, useful for template morphing
00357 
00358   double shiftBy_; // option to increase/decrease within uncertainties the jet energy resolution used for smearing 
00359 
00360   StringCutObjectSelector<T>* skipJetSelection_; // jets passing this cut are **not** smeared 
00361   double skipRawJetPtThreshold_;  // jets with transverse momenta below this value (either on "raw" or "corrected" level) 
00362   double skipCorrJetPtThreshold_; // are **not** smeared 
00363 
00364   int verbosity_; // flag to enabled/disable debug output
00365 };
00366 
00367 #endif