CMS 3D CMS Logo

FFTJetCorrectionProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: RecoJets/FFTJetProducers
4 // Class: FFTJetCorrectionProducer
5 //
13 //
14 // Original Author: Igor Volobouev
15 // Created: Mon Aug 6 11:03:38 CDT 2012
16 //
17 //
18 
19 
20 // system include files
21 #include <iostream>
22 #include <memory>
23 #include <cfloat>
24 #include <cmath>
25 
26 // user include files
29 
32 
34 
38 
39 
40 #define PILEUP_CALCULATION_MASK 0x200
41 #define PILEUP_SUBTRACTION_MASK_4VEC 0x400
42 #define PILEUP_SUBTRACTION_MASK_PT 0x800
43 #define PILEUP_SUBTRACTION_MASK_ANY (PILEUP_SUBTRACTION_MASK_4VEC | \
44  PILEUP_SUBTRACTION_MASK_PT)
45 
46 using namespace fftjetcms;
47 
48 //
49 // A generic switch statement based on jet type
50 //
51 #define jet_type_switch(method, arg1, arg2) do {\
52  switch (jetType)\
53  {\
54  case CALOJET:\
55  method <reco::CaloJet> ( arg1 , arg2 );\
56  break;\
57  case PFJET:\
58  method <reco::PFJet> ( arg1 , arg2 );\
59  break;\
60  case GENJET:\
61  method <reco::GenJet> ( arg1 , arg2 );\
62  break;\
63  case TRACKJET:\
64  method <reco::TrackJet> ( arg1 , arg2 );\
65  break;\
66  case BASICJET:\
67  method <reco::BasicJet> ( arg1 , arg2 );\
68  break;\
69  case JPTJET:\
70  method <reco::JPTJet> ( arg1 , arg2 );\
71  break;\
72  default:\
73  assert(!"ERROR in FFTJetCorrectionProducer : invalid jet type."\
74  " This is a bug. Please report.");\
75  }\
76 } while(0);
77 
78 
79 namespace {
80  struct LocalSortByPt
81  {
82  template<class Jet>
83  inline bool operator()(const Jet& l, const Jet& r) const
84  {return l.pt() > r.pt();}
85  };
86 }
87 
88 
89 //
90 // class declaration
91 //
93 {
94 public:
96  ~FFTJetCorrectionProducer() override;
97 
98 private:
99 
100  void beginJob() override ;
101  void produce(edm::Event&, const edm::EventSetup&) override;
102  void endJob() override ;
103 
104  template <typename Jet>
105  void makeProduces(const std::string& alias, const std::string& tag);
106 
107  template <typename Jet>
108  void applyCorrections(edm::Event& iEvent, const edm::EventSetup& iSetup);
109 
110  template <typename Jet>
111  void performPileupSubtraction(Jet&);
112 
113  // Label for the input collection
115 
116  // Label for the output objects
118 
119  // Jet type to process
121 
122  // The names of the jet correction records
123  const std::vector<std::string> records;
124 
125  // Are we going to create the uncertainty collection?
126  const bool writeUncertainties;
127 
128  // What to do about pileup subtraction
129  const bool subtractPileup;
131 
132  // Print some info about jets
133  const bool verbose;
134 
135  // Space for level masks
136  std::vector<int> sequenceMasks;
137 
138  // Event counter
139  unsigned long eventCount;
140 
141  // tokens for data access
143 };
144 
145 
146 template <typename T>
148  const std::string& alias, const std::string& tag)
149 {
150  produces<std::vector<reco::FFTAnyJet<T> > >(tag).setBranchAlias(alias);
151 }
152 
153 
154 template <typename Jet>
156 {
157  using reco::FFTJet;
158 
159  FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(jet.getFFTSpecific()));
160  const math::XYZTLorentzVector& new4vec = adjustForPileup(
161  fftJet.f_vec(), fftJet.f_pileup(), subtractPileupAs4Vec);
162  fftJet.setFourVec(new4vec);
163  int status = fftJet.f_status();
166  else
167  status |= PILEUP_SUBTRACTION_MASK_PT;
168  fftJet.setStatus(status);
169  jet.setP4(new4vec);
170 }
171 
172 
173 template <typename Jet>
175  const edm::EventSetup& iSetup)
176 {
177  using reco::FFTJet;
178  typedef reco::FFTAnyJet<Jet> MyJet;
179  typedef std::vector<MyJet> MyCollection;
180  typedef typename FFTJetCorrectorSequenceTypemap<MyJet>::loader Loader;
181  typedef typename Loader::data_type CorrectorSequence;
182  typedef typename CorrectorSequence::result_type CorrectorResult;
183 
184  // Load the jet corrector sequences
185  const unsigned nRecords = records.size();
186  std::vector<edm::ESHandle<CorrectorSequence> > handles(nRecords);
187  for (unsigned irec=0; irec<nRecords; ++irec)
188  Loader::instance().load(iSetup, records[irec], handles[irec]);
189 
190  // Figure out which correction levels we are applying
191  // and create masks which will indicate this
192  sequenceMasks.clear();
193  sequenceMasks.reserve(nRecords);
194 
195  int totalMask = 0;
196  for (unsigned irec=0; irec<nRecords; ++irec)
197  {
198  int levelMask = 0;
199  const unsigned nLevels = handles[irec]->nLevels();
200  for (unsigned i=0; i<nLevels; ++i)
201  {
202  const unsigned lev = (*handles[irec])[i].level();
203 
204  // Not tracking "level 0" corrections in the status word.
205  // Level 0 is basically reserved for uncertainty calculations.
206  if (lev)
207  {
208  const int mask = (1 << lev);
209  if (totalMask & mask)
210  throw cms::Exception("FFTJetBadConfig")
211  << "Error in FFTJetCorrectionProducer::applyCorrections:"
212  << " jet correction at level " << lev
213  << " is applied more than once\n";
214  totalMask |= mask;
215  levelMask |= mask;
216  }
217  }
218  sequenceMasks.push_back(levelMask << 12);
219  }
220  totalMask = (totalMask << 12);
221 
222  // Is this data or MC?
223  const bool isMC = !iEvent.isRealData();
224 
225  // Load the jet collection
227  iEvent.getByToken(input_jets_token_, jets);
228 
229  // Create the output collection
230  const unsigned nJets = jets->size();
231  auto coll = std::make_unique<MyCollection>();
232  coll->reserve(nJets);
233 
234  // Cycle over jets and apply the corrector sequences
235  bool sorted = true;
236  double previousPt = DBL_MAX;
237  for (unsigned ijet=0; ijet<nJets; ++ijet)
238  {
239  const MyJet& j((*jets)[ijet]);
240 
241  // Check that this jet has not been corrected yet
242  const int initialStatus = j.getFFTSpecific().f_status();
243  if (initialStatus & totalMask)
244  throw cms::Exception("FFTJetBadConfig")
245  << "Error in FFTJetCorrectionProducer::applyCorrections: "
246  << "this jet collection is already corrected for some or all "
247  << "of the specified levels\n";
248 
249  MyJet corJ(j);
250 
251  if (verbose)
252  {
253  const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
254  std::cout << "++++ Evt " << eventCount << " jet " << ijet
255  << ": pt = " << corJ.pt()
256  << ", eta = " << fj.f_vec().eta()
257  << ", R = " << fj.f_recoScale()
258  << ", s = 0x" << std::hex << fj.f_status() << std::dec
259  << std::endl;
260  }
261 
262  // Check if we need to subtract pileup first.
263  // Pileup subtraction is not part of the corrector sequence
264  // itself because 4-vector subtraction does not map well
265  // into multiplication of 4-vectors by a scale factor.
266  if (subtractPileup)
267  {
268  if (initialStatus & PILEUP_SUBTRACTION_MASK_ANY)
269  throw cms::Exception("FFTJetBadConfig")
270  << "Error in FFTJetCorrectionProducer::applyCorrections: "
271  << "this jet collection is already pileup-subtracted\n";
272  if (!(initialStatus & PILEUP_CALCULATION_MASK))
273  throw cms::Exception("FFTJetBadConfig")
274  << "Error in FFTJetCorrectionProducer::applyCorrections: "
275  << "pileup was not calculated for this jet collection\n";
276  performPileupSubtraction(corJ);
277 
278  if (verbose)
279  {
280  const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
281  std::cout << " Pileup subtract"
282  << ": pt = " << corJ.pt()
283  << ", eta = " << fj.f_vec().eta()
284  << ", R = " << fj.f_recoScale()
285  << ", s = 0x" << std::hex << fj.f_status() << std::dec
286  << std::endl;
287  }
288  }
289 
290  // Apply all jet correction sequences
291  double sigmaSquared = 0.0;
292  for (unsigned irec=0; irec<nRecords; ++irec)
293  {
294  const CorrectorResult& corr = handles[irec]->correct(corJ, isMC);
295 
296  // Update the 4-vector
297  FFTJet<float>& fftJet(const_cast<FFTJet<float>&>(corJ.getFFTSpecific()));
298  corJ.setP4(corr.vec());
299  fftJet.setFourVec(corr.vec());
300 
301  // Update the jet correction status
302  fftJet.setStatus(fftJet.f_status() | sequenceMasks[irec]);
303 
304  // Update the (systematic) uncertainty
305  const double s = corr.sigma();
306  sigmaSquared += s*s;
307  }
308 
309  // There is no place for uncertainty in the jet structure.
310  // However, there is the unused pileup field (FFTJet maintains
311  // the pileup separately as a 4-vector). Use this unused field
312  // to store the uncertainty. This hack is needed because
313  // subsequent sequence sorting by Pt can change the jet ordering.
314  if (writeUncertainties)
315  corJ.setPileup(sqrt(sigmaSquared));
316 
317  coll->push_back(corJ);
318 
319  // Check whether the sequence remains sorted by pt
320  const double pt = corJ.pt();
321  if (pt > previousPt)
322  sorted = false;
323  previousPt = pt;
324 
325  if (verbose)
326  {
327  const reco::FFTJet<float>& fj = corJ.getFFTSpecific();
328  std::cout << " Fully corrected"
329  << ": pt = " << corJ.pt()
330  << ", eta = " << fj.f_vec().eta()
331  << ", R = " << fj.f_recoScale()
332  << ", s = 0x" << std::hex << fj.f_status() << std::dec
333  << std::endl;
334  }
335  }
336 
337  if (!sorted)
338  std::sort(coll->begin(), coll->end(), LocalSortByPt());
339 
340  // Create the uncertainty sequence
341  if (writeUncertainties)
342  {
343  auto unc = std::make_unique<std::vector<float>>();
344  unc->reserve(nJets);
345  for (unsigned ijet=0; ijet<nJets; ++ijet)
346  {
347  MyJet& j((*coll)[ijet]);
348  unc->push_back(j.pileup());
349  j.setPileup(0.f);
350  }
351  iEvent.put(std::move(unc), outputLabel);
352  }
353 
354  iEvent.put(std::move(coll), outputLabel);
355  ++eventCount;
356 }
357 
358 
359 //
360 // constructors and destructor
361 //
363  : inputLabel(ps.getParameter<edm::InputTag>("src")),
364  outputLabel(ps.getParameter<std::string>("outputLabel")),
365  jetType(parseJetType(ps.getParameter<std::string>("jetType"))),
366  records(ps.getParameter<std::vector<std::string> >("records")),
367  writeUncertainties(ps.getParameter<bool>("writeUncertainties")),
368  subtractPileup(ps.getParameter<bool>("subtractPileup")),
369  subtractPileupAs4Vec(ps.getParameter<bool>("subtractPileupAs4Vec")),
370  verbose(ps.getUntrackedParameter<bool>("verbose", false)),
371  eventCount(0UL)
372 {
374  "alias", outputLabel));
376 
377  if (writeUncertainties)
378  produces<std::vector<float> >(outputLabel).setBranchAlias(alias);
379 
380  input_jets_token_ = consumes<std::vector<reco::FFTAnyJet<reco::Jet> > >(inputLabel);
381 }
382 
383 
385 {
386 }
387 
388 
389 // ------------ method called to produce the data ------------
391  const edm::EventSetup& iSetup)
392 {
393  jet_type_switch(applyCorrections, iEvent, iSetup);
394 }
395 
396 // ------------ method called once each job just before starting event loop ------------
398 {
399 }
400 
401 // ------------ method called once each job just after ending the event loop ------------
403 {
404 }
405 
406 //define this as a plug-in
T getUntrackedParameter(std::string const &, T const &) const
edm::EDGetTokenT< std::vector< reco::FFTAnyJet< reco::Jet > > > input_jets_token_
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
static PFTauRenderPlugin instance
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
void push_back(key_type i, value_type const &j, label_type const &flav="")
int f_status() const
Definition: FFTJet.h:77
#define jet_type_switch(method, arg1, arg2)
bool isRealData() const
Definition: EventBase.h:62
void produce(edm::Event &, const edm::EventSetup &) override
void makeProduces(const std::string &alias, const std::string &tag)
void beginJob()
Definition: Breakpoints.cc:14
void applyCorrections(edm::Event &iEvent, const edm::EventSetup &iSetup)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int iEvent
Definition: GenABIO.cc:224
#define PILEUP_SUBTRACTION_MASK_PT
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Real f_recoScale() const
Definition: FFTJet.h:73
Definition: Jet.py:1
T sqrt(T t)
Definition: SSEVec.h:18
vector< PseudoJet > jets
Implements inheritance relationships for FFTJet jets.
Definition: FFTAnyJet.h:16
#define PILEUP_SUBTRACTION_MASK_4VEC
const math::XYZTLorentzVector & f_vec() const
Definition: FFTJet.h:62
double f[11][100]
FFTJetCorrectionProducer(const edm::ParameterSet &)
JetType parseJetType(const std::string &name)
Definition: JetType.cc:5
JetCorrectorParameters corr
Definition: classes.h:5
#define PILEUP_SUBTRACTION_MASK_ANY
Storage class for jets reconstructed by FFTJet package.
Definition: FFTJet.h:19
JetCorrectorParametersCollection coll
Definition: classes.h:10
HLT enums.
#define PILEUP_CALCULATION_MASK
const std::vector< std::string > records
math::XYZTLorentzVector adjustForPileup(const math::XYZTLorentzVector &jet, const math::XYZTLorentzVector &pileup, bool subtractPileupAs4Vec)
def move(src, dest)
Definition: eostools.py:511