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
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
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
00117
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
00154
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
00179
00180 shiftBy_ = ( cfg.exists("shiftBy") ) ? cfg.getParameter<double>("shiftBy") : 0.;
00181
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
00264
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
00277
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
00286
00287
00288
00289
00290
00291 smearedJetEn = jet.energy()*(1. + rnd_.Gaus(0., sigmaEn)/TMath::Max(rawJetP4.E(), corrJetP4.E()));
00292 }
00293 }
00294
00295
00296 const double minJetEn = 1.e-2;
00297 if ( smearedJetEn < minJetEn ) smearedJetEn = minJetEn;
00298
00299
00300
00301
00302
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
00326 evt.put(smearedJets);
00327 }
00328
00329 std::string moduleLabel_;
00330
00331 SmearedJetProducer_namespace::GenJetMatcherT<T> genJetMatcher_;
00332
00333
00334
00335
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_;
00345 double jetCorrEtaMax_;
00346
00347
00348
00349
00350 Textractor jetCorrExtractor_;
00351
00352 double sigmaMaxGenJetMatch_;
00353
00354
00355
00356 double smearBy_;
00357
00358 double shiftBy_;
00359
00360 StringCutObjectSelector<T>* skipJetSelection_;
00361 double skipRawJetPtThreshold_;
00362 double skipCorrJetPtThreshold_;
00363
00364 int verbosity_;
00365 };
00366
00367 #endif