CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SmearedJetProducerT.h
Go to the documentation of this file.
1 #ifndef PhysicsTools_PatUtils_SmearedJetProducerT_h
2 #define PhysicsTools_PatUtils_SmearedJetProducerT_h
3 
24 
27 
31 
35 
39 
40 #include <TFile.h>
41 #include <TFormula.h>
42 #include <TH2.h>
43 #include <TRandom3.h>
44 #include <TString.h>
45 
46 #include <iostream>
47 #include <iomanip>
48 #include <type_traits>
49 
50 namespace SmearedJetProducer_namespace
51 {
52  template <typename T>
54  {
55  public:
56 
58  : srcGenJetsToken_(iC.consumes<reco::GenJetCollection>(cfg.getParameter<edm::InputTag>("srcGenJets"))),
60  {
61  TString dRmaxGenJetMatch_formula = cfg.getParameter<std::string>("dRmaxGenJetMatch").data();
62  dRmaxGenJetMatch_formula.ReplaceAll("genJetPt", "x");
63  dRmaxGenJetMatch_ = new TFormula("dRmaxGenJetMatch", dRmaxGenJetMatch_formula.Data());
64  }
66  {
67  delete dRmaxGenJetMatch_;
68  }
69 
70  const reco::GenJet* operator()(const T& jet, edm::Event* evt = 0) const
71  {
72  assert(evt);
73 
75  evt->getByToken(srcGenJetsToken_, genJets);
76 
77  const reco::GenJet* retVal = 0;
78 
79  double dRbestMatch = 1.e+6;
80  for ( reco::GenJetCollection::const_iterator genJet = genJets->begin();
81  genJet != genJets->end(); ++genJet ) {
82  double dRmax = dRmaxGenJetMatch_->Eval(genJet->pt());
83  //std::cout << "genJetPt = " << genJet->pt() << ": dRmax = " << dRmax << std::endl;
84  double dR = deltaR(jet.p4(), genJet->p4());
85  if ( dR < dRbestMatch && dR < dRmax ) {
86  retVal = &(*genJet);
87  dRbestMatch = dR;
88  }
89  }
90 
91  return retVal;
92  }
93 
94  private:
95 
96 //--- configuration parameter
99 
100  TFormula* dRmaxGenJetMatch_;
101  };
102 
103  template <typename T>
105  {
106  public:
107 
110 
111  double operator()(const T&) const
112  {
113  throw cms::Exception("SmearedJetProducer::produce")
114  << " Jets of type other than PF not supported yet !!\n";
115  }
116  };
117 
118  template <typename T>
119  class RawJetExtractorT // this template is neccessary to support pat::Jets
120  // (because pat::Jet->p4() returns the JES corrected, not the raw, jet momentum)
121  {
122  public:
123 
126  {
127  return jet.p4();
128  }
129  };
130 
131  template <>
132  class RawJetExtractorT<pat::Jet>
133  {
134  public:
135 
138  {
139  if ( jet.jecSetsAvailable() ) return jet.correctedP4("Uncorrected");
140  else return jet.p4();
141  }
142  };
143 }
144 
145 template <typename T, typename Textractor>
147 {
148  typedef std::vector<T> JetCollection;
149 
150  public:
151 
153  : moduleLabel_(cfg.getParameter<std::string>("@module_label")),
155  jetResolutionExtractor_(cfg.getParameter<edm::ParameterSet>("jetResolutions")),
156  jetCorrLabel_(""),
158  {
159  //std::cout << "<SmearedJetProducer::SmearedJetProducer>:" << std::endl;
160  //std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
161 
162  src_ = cfg.getParameter<edm::InputTag>("src");
163  srcToken_ = consumes<JetCollection>(src_);
164 
166  std::string lutName = cfg.getParameter<std::string>("lutName");
167  if (inputFileName.location() == edm::FileInPath::Unknown)
168  throw cms::Exception("JetMETsmearInputProducer")
169  << " Failed to find File = " << inputFileName << " !!\n";
170 
171  inputFile_ = new TFile(inputFileName.fullPath().data());
172  lut_ = dynamic_cast<TH2*>(inputFile_->Get(lutName.data()));
173  if ( !lut_ )
174  throw cms::Exception("SmearedJetProducer")
175  << " Failed to load LUT = " << lutName.data() << " from file = " << inputFileName.fullPath().data() << " !!\n";
176 
177  if ( cfg.exists("jetCorrLabel") ) {
178  jetCorrLabel_ = cfg.getParameter<edm::InputTag>("jetCorrLabel");
179  jetCorrToken_ = consumes<reco::JetCorrector>(jetCorrLabel_);
180  }
181 
182  jetCorrEtaMax_ = ( cfg.exists("jetCorrEtaMax") ) ?
183  cfg.getParameter<double>("jetCorrEtaMax") : 9.9;
184 
185  sigmaMaxGenJetMatch_ = cfg.getParameter<double>("sigmaMaxGenJetMatch");
186 
187  smearBy_ = ( cfg.exists("smearBy") ) ? cfg.getParameter<double>("smearBy") : 1.0;
188  //std::cout << "smearBy = " << smearBy_ << std::endl;
189 
190  shiftBy_ = ( cfg.exists("shiftBy") ) ? cfg.getParameter<double>("shiftBy") : 0.;
191  //std::cout << "shiftBy = " << shiftBy_ << std::endl;
192 
193  if ( cfg.exists("skipJetSelection") ) {
194  std::string skipJetSelection_string = cfg.getParameter<std::string>("skipJetSelection");
195  skipJetSelection_ = new StringCutObjectSelector<T>(skipJetSelection_string);
196  }
197 
198  skipRawJetPtThreshold_ = ( cfg.exists("skipRawJetPtThreshold") ) ?
199  cfg.getParameter<double>("skipRawJetPtThreshold") : 1.e-2;
200  skipCorrJetPtThreshold_ = ( cfg.exists("skipCorrJetPtThreshold") ) ?
201  cfg.getParameter<double>("skipCorrJetPtThreshold") : 1.e-2;
202 
203  verbosity_ = ( cfg.exists("verbosity") ) ?
204  cfg.getParameter<int>("verbosity") : 0;
205 
206  produces<JetCollection>();
207  }
209  {
210  delete skipJetSelection_;
211  delete inputFile_;
212  delete lut_;
213  }
214 
215  private:
216 
217  virtual void produce(edm::Event& evt, const edm::EventSetup& es)
218  {
219  if ( verbosity_ ) {
220  std::cout << "<SmearedJetProducerT::produce>:" << std::endl;
221  std::cout << " moduleLabel = " << moduleLabel_ << std::endl;
222  std::cout << " src = " << src_.label() << std::endl;
223  }
224 
225  std::auto_ptr<JetCollection> smearedJets(new JetCollection);
226 
228  evt.getByToken(srcToken_, jets);
229 
230  int numJets = jets->size();
231  for ( int jetIndex = 0; jetIndex < numJets; ++jetIndex ) {
232  const T& jet = jets->at(jetIndex);
233 
234  static const SmearedJetProducer_namespace::RawJetExtractorT<T> rawJetExtractor;
235  reco::Candidate::LorentzVector rawJetP4 = rawJetExtractor(jet);
236  if ( verbosity_ ) {
237  std::cout << "rawJet: Pt = " << rawJetP4.pt() << ", eta = " << rawJetP4.eta() << ", phi = " << rawJetP4.phi() << std::endl;
238  }
239 
240  reco::Candidate::LorentzVector corrJetP4 = jet.p4();
241  if ( !jetCorrLabel_.label().empty() ) {
243  evt.getByToken(jetCorrToken_, jetCorr);
244  corrJetP4 =
247  jetCorrExtractor_(jet, jetCorr.product(), jetCorrEtaMax_, &rawJetP4);
248  }
249  if ( verbosity_ ) {
250  std::cout << "corrJet: Pt = " << corrJetP4.pt() << ", eta = " << corrJetP4.eta() << ", phi = " << corrJetP4.phi() << std::endl;
251  }
252 
253  double smearFactor = 1.;
254  double x = std::abs(corrJetP4.eta());
255  double y = corrJetP4.pt();
256  if ( x > lut_->GetXaxis()->GetXmin() && x < lut_->GetXaxis()->GetXmax() &&
257  y > lut_->GetYaxis()->GetXmin() && y < lut_->GetYaxis()->GetXmax() ) {
258  int binIndex = lut_->FindBin(x, y);
259 
260  if ( smearBy_ > 0. ) smearFactor += smearBy_*(lut_->GetBinContent(binIndex) - 1.);
261  double smearFactorErr = lut_->GetBinError(binIndex);
262  if ( verbosity_ ) std::cout << "smearFactor = " << smearFactor << " +/- " << smearFactorErr << std::endl;
263 
264  if ( shiftBy_ != 0. ) {
265  smearFactor += (shiftBy_*smearFactorErr);
266  if ( verbosity_ ) std::cout << "smearFactor(shifted) = " << smearFactor << std::endl;
267  }
268  }
269 
270  double smearedJetEn = jet.energy();
271  double sigmaEn = jetResolutionExtractor_(jet)*sqrt(smearFactor*smearFactor - 1.);
272  const reco::GenJet* genJet = genJetMatcher_(jet, &evt);
273  bool isGenMatched = false;
274  if ( genJet ) {
275  if ( verbosity_ ) {
276  std::cout << "genJet: Pt = " << genJet->pt() << ", eta = " << genJet->eta() << ", phi = " << genJet->phi() << std::endl;
277  }
278  double dEn = corrJetP4.E() - genJet->energy();
279  if ( dEn < (sigmaMaxGenJetMatch_*sigmaEn) ) {
280 //--- case 1: reconstructed jet matched to generator level jet,
281 // smear difference between reconstructed and "true" jet energy
282 
283  if ( verbosity_ ) {
284  std::cout << " successfully matched to genJet" << std::endl;
285  std::cout << "corrJetEn = " << corrJetP4.E() << ", genJetEn = " << genJet->energy() << " --> dEn = " << dEn << std::endl;
286  }
287 
288  smearedJetEn = jet.energy()*(1. + (smearFactor - 1.)*dEn/std::max( rawJetP4.E(), corrJetP4.E()));
289  isGenMatched = true;
290  }
291  }
292  if ( !isGenMatched ) {
293 //--- case 2: reconstructed jet **not** matched to generator level jet,
294 // smear jet energy using MC resolution functions implemented in PFMEt significance algorithm (CMS AN-10/400)
295 
296  if ( verbosity_ ) {
297  std::cout << " not matched to genJet" << std::endl;
298  std::cout << "corrJetEn = " << corrJetP4.E() << ", sigmaEn = " << sigmaEn << std::endl;
299  }
300 
301  if ( smearFactor > 1. ) {
302  // CV: MC resolution already accounted for in reconstructed jet,
303  // add additional Gaussian smearing of width = sqrt(smearFactor^2 - 1)
304  // to account for Data/MC **difference** in jet resolutions.
305  // Take maximum(rawJetEn, corrJetEn) to avoid pathological cases
306  // (e.g. corrJetEn << rawJetEn, due to L1Fastjet corrections)
307 
308  smearedJetEn = jet.energy()*(1. + rnd_.Gaus(0., sigmaEn)/std::max( rawJetP4.E(), corrJetP4.E()) );
309  }
310  }
311 
312  // CV: keep minimum jet energy, in order not to loose direction information
313  const double minJetEn = 1.e-2;
314  if ( smearedJetEn < minJetEn ) smearedJetEn = minJetEn;
315 
316  // CV: skip smearing in case either "raw" or "corrected" jet energy is very low
317  // or jet passes selection configurable via python
318  // (allows for protection against "pathological cases",
319  // cf. PhysicsTools/PatUtils/python/tools/metUncertaintyTools.py)
320  reco::Candidate::LorentzVector smearedJetP4 = jet.p4();
321  if ( !((skipJetSelection_ && (*skipJetSelection_)(jet)) ||
322  rawJetP4.pt() < skipRawJetPtThreshold_ ||
323  corrJetP4.pt() < skipCorrJetPtThreshold_ ) ) {
324  if ( verbosity_ ) {
325  std::cout << " smearing jetP4 by factor = " << (smearedJetEn/jet.energy()) << " --> smearedJetEn = " << smearedJetEn << std::endl;
326  }
327  smearedJetP4 *= (smearedJetEn/jet.energy());
328  }
329 
330  if ( verbosity_ ) {
331  std::cout << "smearedJet: Pt = " << smearedJetP4.pt() << ", eta = " << smearedJetP4.eta() << ", phi = " << smearedJetP4.phi() << std::endl;
332  std::cout << " dPt = " << (smearedJetP4.pt() - jet.pt())
333  << " (Px = " << (smearedJetP4.px() - jet.px()) << ", Py = " << (smearedJetP4.py() - jet.py()) << ")" << std::endl;
334  }
335 
336  T smearedJet = (jet);
337  smearedJet.setP4(smearedJetP4);
338 
339  smearedJets->push_back(smearedJet);
340  }
341 
342 //--- add collection of "smeared" jets to the event
343  evt.put(smearedJets);
344  }
345 
347 
349 
350 //--- configuration parameters
351 
352  // collection of pat::Jets (with L2L3/L2L3Residual corrections applied)
355 
356  TFile* inputFile_;
357  TH2* lut_;
358 
360  TRandom3 rnd_;
361 
363  edm::EDGetTokenT<reco::JetCorrector> jetCorrToken_; // e.g. 'ak5CaloJetL1FastL2L3' (MC) / 'ak5CaloJetL1FastL2L3Residual' (Data)
364  double jetCorrEtaMax_; // do not use JEC factors for |eta| above this threshold (recommended default = 4.7),
365  // in order to work around problem with CMSSW_4_2_x JEC factors at high eta,
366  // reported in
367  // https://hypernews.cern.ch/HyperNews/CMS/get/jes/270.html
368  // https://hypernews.cern.ch/HyperNews/CMS/get/JetMET/1259/1.html
369  Textractor jetCorrExtractor_;
370 
371  double sigmaMaxGenJetMatch_; // maximum difference between energy of reconstructed jet and matched generator level jet
372  // (if the difference between reconstructed and generated jet energy exceeds this threshold,
373  // the jet is considered to have substantial pile-up contributions are is considered to be unmatched)
374 
375  double smearBy_; // option to "smear" jet energy by N standard-deviations, useful for template morphing
376 
377  double shiftBy_; // option to increase/decrease within uncertainties the jet energy resolution used for smearing
378 
379  StringCutObjectSelector<T>* skipJetSelection_; // jets passing this cut are **not** smeared
380  double skipRawJetPtThreshold_; // jets with transverse momenta below this value (either on "raw" or "corrected" level)
381  double skipCorrJetPtThreshold_; // are **not** smeared
382 
383  int verbosity_; // flag to enabled/disable debug output
384 };
385 
386 #endif
T getParameter(std::string const &) const
reco::Candidate::LorentzVector operator()(const pat::Jet &jet) const
edm::EDGetTokenT< JetCollection > srcToken_
StringCutObjectSelector< T > * skipJetSelection_
bool jecSetsAvailable() const
Definition: Jet.h:122
virtual float pt() const
transverse momentum
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
virtual float phi() const
momentum azimuthal angle
std::vector< GenJet > GenJetCollection
collection of GenJet objects
bool exists(std::string const &parameterName) const
checks if a parameter exists
const LorentzVector & correctedP4(const std::string &level, const std::string &flavor="none", const std::string &set="") const
Definition: Jet.h:149
SmearedJetProducer_namespace::GenJetMatcherT< T > genJetMatcher_
SmearedJetProducerT(const edm::ParameterSet &cfg)
virtual double energy() const
energy
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:113
virtual float eta() const
momentum pseudorapidity
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T sqrt(T t)
Definition: SSEVec.h:48
vector< PseudoJet > jets
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Jets made from MC generator particles.
Definition: GenJet.h:24
edm::EDGetTokenT< reco::JetCorrector > jetCorrToken_
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
T const * product() const
Definition: Handle.h:81
GenJetMatcherT(const edm::ParameterSet &cfg, edm::ConsumesCollector &&iC)
SmearedJetProducer_namespace::JetResolutionExtractorT< T > jetResolutionExtractor_
Analysis-level calorimeter jet class.
Definition: Jet.h:76
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:41
std::string const & label() const
Definition: InputTag.h:42
const reco::GenJet * operator()(const T &jet, edm::Event *evt=0) const
edm::EDGetTokenT< reco::GenJetCollection > srcGenJetsToken_
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
reco::Candidate::LorentzVector operator()(const T &jet) const
tuple cout
Definition: gather_cfg.py:121
Definition: DDAxes.h:10
long double T
std::vector< T > JetCollection
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
virtual void produce(edm::Event &evt, const edm::EventSetup &es)