CMS 3D CMS Logo

L1JetRecoTreeProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: L1Trigger/L1TNtuples
4 // Class: L1JetRecoTreeProducer
5 //
13 // system include files
14 #include <memory>
15 
16 // framework
25 
26 // cond formats
30 
31 // data formats
41 
42 // ROOT output stuff
45 #include "TH1.h"
46 #include "TTree.h"
47 #include "TF1.h"
48 #include <TVector2.h>
49 
50 //local data formats
53 
54 //
55 // class declaration
56 //
57 
58 class L1JetRecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
59 public:
61  ~L1JetRecoTreeProducer() override;
62 
63 private:
64  void beginJob(void) override;
65  void analyze(const edm::Event&, const edm::EventSetup&) override;
66  void endJob() override;
67 
74 
77 
78  bool pfJetID(const reco::PFJet& jet);
79  bool caloJetID(const reco::CaloJet& jet);
80 
81 public:
84 
85 private:
86  // tree
87  TTree* tree_;
88 
89  // EDM input tags
95 
99 
101 
102  // debug stuff
105  double jetetaMax_;
106  unsigned int maxCl_;
107  unsigned int maxJet_;
108  unsigned int maxVtx_;
109  unsigned int maxTrk_;
110 
118 
120 };
121 
123  : pfJetsMissing_(false),
124  pfJetCorrMissing_(false),
125  caloJetCorrMissing_(false),
126  caloJetsMissing_(false),
127  caloJetIDMissing_(false),
128  pfMetMissing_(false),
129  caloMetMissing_(false),
130  caloMetBEMissing_(false),
131  muonsMissing_(false) {
132  caloJetToken_ =
133  consumes<reco::CaloJetCollection>(iConfig.getUntrackedParameter("caloJetToken", edm::InputTag("ak4CaloJets")));
134  pfJetToken_ =
135  consumes<reco::PFJetCollection>(iConfig.getUntrackedParameter("pfJetToken", edm::InputTag("ak4PFJetsCHS")));
137  consumes<edm::ValueMap<reco::JetID> >(iConfig.getUntrackedParameter("caloJetIDToken", edm::InputTag("ak4JetID")));
138  pfJECToken_ = consumes<reco::JetCorrector>(
139  iConfig.getUntrackedParameter<edm::InputTag>("pfJECToken", edm::InputTag("ak4PFCHSL1FastL2L3ResidualCorrector")));
140  caloJECToken_ = consumes<reco::JetCorrector>(iConfig.getUntrackedParameter<edm::InputTag>(
141  "caloJECToken", edm::InputTag("ak4CaloL1FastL2L3ResidualCorrector")));
142 
143  pfMetToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("pfMetToken", edm::InputTag("pfMetT1")));
144  caloMetToken_ =
145  consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetToken", edm::InputTag("caloMet")));
147  consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetBEToken", edm::InputTag("caloMetBE")));
148 
149  muonToken_ = consumes<reco::MuonCollection>(iConfig.getUntrackedParameter("muonToken", edm::InputTag("muons")));
150 
151  usesResource(TFileService::kSharedResource);
152 
153  jetptThreshold_ = iConfig.getParameter<double>("jetptThreshold");
154  jetetaMax_ = iConfig.getParameter<double>("jetetaMax");
155  maxJet_ = iConfig.getParameter<unsigned int>("maxJet");
156 
159 
160  // set up output
162  tree_ = fs_->make<TTree>("JetRecoTree", "JetRecoTree");
163  tree_->Branch("Jet", "L1Analysis::L1AnalysisRecoJetDataFormat", &jet_data, 32000, 3);
164  tree_->Branch("Sums", "L1Analysis::L1AnalysisRecoMetDataFormat", &met_data, 32000, 3);
165 }
166 
168  // do anything here that needs to be done at desctruction time
169  // (e.g. close files, deallocate resources etc.)
170 }
171 
172 //
173 // member functions
174 //
175 
176 // ------------ method called to for each event ------------
178  jet_data->Reset();
179  met_data->Reset();
180 
181  // get jets
183  iEvent.getByToken(pfJetToken_, pfJets);
184 
185  // get calo jets
187  iEvent.getByToken(caloJetToken_, caloJets);
188 
189  //get sums
191  iEvent.getByToken(pfMetToken_, pfMet);
192 
193  // get jet ID
195  iEvent.getByToken(caloJetIDToken_, jetsID);
196 
198  iEvent.getByToken(pfJECToken_, pfJetCorr);
199 
201  iEvent.getByToken(caloJECToken_, caloJetCorr);
202 
204  iEvent.getByToken(caloMetToken_, caloMet);
205 
207  iEvent.getByToken(caloMetBEToken_, caloMetBE);
208 
209  // get muons
211  iEvent.getByToken(muonToken_, muons);
212 
213  if (pfJets.isValid()) {
214  jet_data->nJets = 0;
215 
216  doPFJets(pfJets);
217 
218  } else {
219  if (!pfJetsMissing_) {
220  edm::LogWarning("MissingProduct") << "PFJets not found. Branch will not be filled" << std::endl;
221  }
222  pfJetsMissing_ = true;
223  }
224 
225  if (pfJetCorr.isValid()) {
226  doPFJetCorr(pfJets, pfJetCorr);
227 
228  } else {
229  if (!pfJetCorrMissing_) {
230  edm::LogWarning("MissingProduct") << "PF Jet Corrector not found. Branch will not be filled" << std::endl;
231  }
232  pfJetCorrMissing_ = true;
233  }
234 
235  if (caloJets.isValid()) {
236  jet_data->nCaloJets = 0;
237 
239 
240  } else {
241  if (!caloJetsMissing_) {
242  edm::LogWarning("MissingProduct") << "Calo Jets not found. Branch will not be filled" << std::endl;
243  }
244  caloJetsMissing_ = true;
245  }
246 
247  if (caloJetCorr.isValid()) {
248  doCaloJetCorr(caloJets, caloJetCorr);
249 
250  } else {
251  if (!caloJetCorrMissing_) {
252  edm::LogWarning("MissingProduct") << "Calo Jet Corrector not found. Branch will not be filled" << std::endl;
253  }
254  caloJetCorrMissing_ = true;
255  }
256 
257  if (!jetsID.isValid()) {
258  if (!caloJetIDMissing_) {
259  edm::LogWarning("MissingProduct") << "Calo Jet ID not found. Branch will not be filled" << std::endl;
260  }
261  caloJetIDMissing_ = true;
262  }
263 
264  if (pfMet.isValid()) {
265  doPFMet(pfMet);
266 
267  if (muons.isValid()) {
269 
270  } else {
271  if (!muonsMissing_) {
272  edm::LogWarning("MissingProduct") << "Muons not found. PFMetNoMu branch will not be filled" << std::endl;
273  }
274  muonsMissing_ = true;
275  }
276  } else {
277  if (!pfMetMissing_) {
278  edm::LogWarning("MissingProduct") << "PFMet not found. Branch will not be filled" << std::endl;
279  }
280  pfMetMissing_ = true;
281  }
282 
283  if (caloMet.isValid()) {
285 
286  } else {
287  if (!caloMetMissing_) {
288  edm::LogWarning("MissingProduct") << "CaloMet not found. Branch will not be filled" << std::endl;
289  }
290  caloMetMissing_ = true;
291  }
292 
293  if (caloMetBE.isValid()) {
295 
296  } else {
297  if (!caloMetBEMissing_) {
298  edm::LogWarning("MissingProduct") << "CaloMetBE not found. Branch will not be filled" << std::endl;
299  }
300  caloMetBEMissing_ = true;
301  }
302 
303  tree_->Fill();
304 }
305 
307  for (auto it = caloJets->begin(); it != caloJets->end() && jet_data->nCaloJets < maxJet_; ++it) {
308  if (!caloJetIDMissing_)
309  if (!caloJetID(*it))
310  continue;
311 
312  jet_data->caloEt.push_back(it->et());
313  jet_data->caloEta.push_back(it->eta());
314  jet_data->caloPhi.push_back(it->phi());
315  jet_data->caloE.push_back(it->energy());
316 
317  jet_data->eEMF.push_back(it->emEnergyFraction());
318  jet_data->eEmEB.push_back(it->emEnergyInEB());
319  jet_data->eEmEE.push_back(it->emEnergyInEE());
320  jet_data->eEmHF.push_back(it->emEnergyInHF());
321  jet_data->eHadHB.push_back(it->hadEnergyInHB());
322  jet_data->eHadHE.push_back(it->hadEnergyInHE());
323  jet_data->eHadHO.push_back(it->hadEnergyInHO());
324  jet_data->eHadHF.push_back(it->hadEnergyInHF());
325  jet_data->eMaxEcalTow.push_back(it->maxEInEmTowers());
326  jet_data->eMaxHcalTow.push_back(it->maxEInHadTowers());
327  jet_data->towerArea.push_back(it->towersArea());
328  jet_data->n60.push_back(it->n60());
329 
330  jet_data->nCaloJets++;
331  }
332 }
333 
335  for (auto it = pfJets->begin(); it != pfJets->end() && jet_data->nJets < maxJet_; ++it) {
336  if (!pfJetID(*it))
337  continue;
338 
339  jet_data->et.push_back(it->et());
340  jet_data->eta.push_back(it->eta());
341  jet_data->phi.push_back(it->phi());
342  jet_data->e.push_back(it->energy());
343 
344  jet_data->chef.push_back(it->chargedHadronEnergyFraction());
345  jet_data->nhef.push_back(it->neutralHadronEnergyFraction());
346  jet_data->pef.push_back(it->photonEnergyFraction());
347  jet_data->eef.push_back(it->electronEnergyFraction());
348  jet_data->mef.push_back(it->muonEnergyFraction());
349  jet_data->hfhef.push_back(it->HFHadronEnergyFraction());
350  jet_data->hfemef.push_back(it->HFEMEnergyFraction());
351  jet_data->chMult.push_back(it->chargedHadronMultiplicity());
352  jet_data->nhMult.push_back(it->neutralHadronMultiplicity());
353  jet_data->phMult.push_back(it->photonMultiplicity());
354  jet_data->elMult.push_back(it->electronMultiplicity());
355  jet_data->muMult.push_back(it->muonMultiplicity());
356  jet_data->hfhMult.push_back(it->HFHadronMultiplicity());
357  jet_data->hfemMult.push_back(it->HFEMMultiplicity());
358 
359  jet_data->cemef.push_back(it->chargedEmEnergyFraction());
360  jet_data->cmef.push_back(it->chargedMuEnergyFraction());
361  jet_data->nemef.push_back(it->neutralEmEnergyFraction());
362  jet_data->cMult.push_back(it->chargedMultiplicity());
363  jet_data->nMult.push_back(it->neutralMultiplicity());
364 
365  jet_data->nJets++;
366  }
367 }
368 
371  float corrFactor = 1.;
372  unsigned int nJets = 0;
373 
374  float mHx = 0;
375  float mHy = 0;
376 
377  met_data->Ht = 0;
378  met_data->mHt = -999.;
379  met_data->mHtPhi = -999.;
380 
381  for (auto it = pfJets->begin(); it != pfJets->end() && nJets < maxJet_; ++it) {
382  if (!pfJetID(*it))
383  continue;
384 
385  corrFactor = pfJetCorr.product()->correction(*it);
386 
387  jet_data->etCorr.push_back(it->et() * corrFactor);
388  jet_data->corrFactor.push_back(corrFactor);
389 
390  nJets++;
391 
392  if (it->pt() * corrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
393  mHx += -1. * it->px() * corrFactor;
394  mHy += -1. * it->py() * corrFactor;
395  met_data->Ht += it->pt() * corrFactor;
396  }
397  }
398 
399  TVector2 tv2 = TVector2(mHx, mHy);
400  met_data->mHt = tv2.Mod();
401  met_data->mHtPhi = tv2.Phi();
402 }
403 
405  edm::Handle<reco::JetCorrector> caloJetCorr) {
406  float caloCorrFactor = 1.;
407  unsigned int nCaloJets = 0;
408 
409  met_data->caloHt = 0;
410 
411  for (auto it = caloJets->begin(); it != caloJets->end() && nCaloJets < maxJet_; ++it) {
412  if (!caloJetIDMissing_)
413  if (!caloJetID(*it))
414  continue;
415 
416  caloCorrFactor = caloJetCorr.product()->correction(*it);
417 
418  jet_data->caloEtCorr.push_back(it->et() * caloCorrFactor);
419  jet_data->caloCorrFactor.push_back(caloCorrFactor);
420 
421  nCaloJets++;
422 
423  if (it->pt() * caloCorrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
424  met_data->caloHt += it->pt() * caloCorrFactor;
425  }
426  }
427 }
428 
430  const reco::PFMETCollection* metCol = pfMet.product();
431  const reco::PFMET theMet = metCol->front();
432 
433  met_data->met = theMet.et();
434  met_data->metPhi = theMet.phi();
435  met_data->sumEt = theMet.sumEt();
436  met_data->metPx = theMet.px();
437  met_data->metPy = theMet.py();
438 }
439 
442  const reco::PFMETCollection* metCol = pfMet.product();
443  const reco::PFMET theMet = metCol->front();
444  reco::PFMET thePFMetNoMu = metCol->front();
445 
446  double pfMetNoMuPx = theMet.px();
447  double pfMetNoMuPy = theMet.py();
448 
449  double muPx(0.), muPy(0.);
450 
451  for (auto it = muons->begin(); it != muons->end(); ++it) {
452  if (it->isPFMuon()) {
453  muPx += it->px();
454  muPy += it->py();
455  }
456  }
457 
458  pfMetNoMuPx += muPx;
459  pfMetNoMuPy += muPy;
460 
461  math::XYZTLorentzVector pfMetNoMuP4(pfMetNoMuPx, pfMetNoMuPy, 0, hypot(pfMetNoMuPx, pfMetNoMuPy));
462 
463  thePFMetNoMu.setP4(pfMetNoMuP4);
464 
465  met_data->pfMetNoMu = thePFMetNoMu.et();
466  met_data->pfMetNoMuPhi = thePFMetNoMu.phi();
467  met_data->pfMetNoMuPx = thePFMetNoMu.px();
468  met_data->pfMetNoMuPy = thePFMetNoMu.py();
469 }
470 
472  const reco::CaloMETCollection* metCol = caloMet.product();
473  const reco::CaloMET theMet = metCol->front();
474 
475  met_data->caloMet = theMet.et();
476  met_data->caloMetPhi = theMet.phi();
477  met_data->caloSumEt = theMet.sumEt();
478 }
479 
481  const reco::CaloMETCollection* metCol = caloMetBE.product();
482  const reco::CaloMET theMet = metCol->front();
483 
484  met_data->caloMetBE = theMet.et();
485  met_data->caloMetPhiBE = theMet.phi();
486  met_data->caloSumEtBE = theMet.sumEt();
487 }
488 
490  bool tmp = true;
491  if (std::abs(jet.eta()) < 2.7) {
492  tmp &= jet.neutralHadronEnergyFraction() < 0.9;
493  tmp &= jet.neutralEmEnergyFraction() < 0.9;
494  tmp &= (jet.chargedMultiplicity() + jet.neutralMultiplicity()) > 1;
495  tmp &= jet.muonEnergyFraction() < 0.8;
496  tmp &= jet.chargedHadronEnergyFraction() > 0.0;
497  tmp &= jet.chargedMultiplicity() > 0;
498  tmp &= jet.chargedEmEnergyFraction() < 0.9;
499  }
500  if (std::abs(jet.eta()) > 2.7 && std::abs(jet.eta()) < 3.0) {
501  tmp &= jet.neutralEmEnergyFraction() > 0.01;
502  tmp &= jet.neutralHadronEnergyFraction() < 0.98;
503  tmp &= jet.neutralMultiplicity() > 2;
504  }
505  if (std::abs(jet.eta()) > 3.0) {
506  tmp &= jet.neutralEmEnergyFraction() < 0.9;
507  tmp &= jet.neutralMultiplicity() > 10;
508  }
509 
510  // our custom selection
511  //tmp &= jet.muonMultiplicity() == 0;
512  //tmp &= jet.electronMultiplicity() == 0;
513 
514  return tmp;
515 }
516 
518  bool tmp = true;
519 
520  return tmp;
521 }
522 
523 // ------------ method called once each job just before starting event loop ------------
525 
526 // ------------ method called once each job just after ending the event loop ------------
528 
529 //define this as a plug-in
static const std::string kSharedResource
Definition: TFileService.h:76
void doPFMet(edm::Handle< reco::PFMETCollection > pfMet)
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
bool pfJetID(const reco::PFJet &jet)
Jets made from CaloTowers.
Definition: CaloJet.h:27
L1Analysis::L1AnalysisRecoMetDataFormat * met_data
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
T const * product() const
Definition: Handle.h:70
void doPFJetCorr(edm::Handle< reco::PFJetCollection > pfJets, edm::Handle< reco::JetCorrector > pfJetCorr)
edm::EDGetTokenT< reco::CaloMETCollection > caloMetToken_
void analyze(const edm::Event &, const edm::EventSetup &) override
void doPFJets(edm::Handle< reco::PFJetCollection > pfJets)
caloMetBE
____________________________________________________________________________||
Definition: CaloMET_cfi.py:19
double sumEt() const
Definition: MET.h:56
edm::EDGetTokenT< reco::JetCorrector > pfJECToken_
double correction(const LorentzVector &fJet) const
get correction using Jet information only
Definition: JetCorrector.h:46
void doCaloMetBE(edm::Handle< reco::CaloMETCollection > caloMetBE)
Jets made from PFObjects.
Definition: PFJet.h:20
T getUntrackedParameter(std::string const &, T const &) const
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
L1Analysis::L1AnalysisRecoJetDataFormat * jet_data
double px() const final
x coordinate of momentum vector
int iEvent
Definition: GenABIO.cc:224
void beginJob(void) override
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::EDGetTokenT< edm::ValueMap< reco::JetID > > caloJetIDToken_
edm::EDGetTokenT< reco::PFJetCollection > pfJetToken_
void doCaloJetCorr(edm::Handle< reco::CaloJetCollection > caloJets, edm::Handle< reco::JetCorrector > caloJetCorr)
double py() const final
y coordinate of momentum vector
void doPFMetNoMu(edm::Handle< reco::PFMETCollection > pfMet, edm::Handle< reco::MuonCollection >)
void doCaloMet(edm::Handle< reco::CaloMETCollection > caloMet)
edm::EDGetTokenT< reco::PFMETCollection > pfMetToken_
std::vector< reco::CaloMET > CaloMETCollection
collection of CaloMET objects
L1JetRecoTreeProducer(const edm::ParameterSet &)
void doCaloJets(edm::Handle< reco::CaloJetCollection > caloJets)
bool caloJetID(const reco::CaloJet &jet)
bool isValid() const
Definition: HandleBase.h:70
edm::EDGetTokenT< reco::MuonCollection > muonToken_
edm::EDGetTokenT< reco::JetCorrector > caloJECToken_
double et() const final
transverse energy
std::vector< reco::PFMET > PFMETCollection
collection of PFMET objects
edm::EDGetTokenT< reco::CaloJetCollection > caloJetToken_
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
Log< level::Warning, false > LogWarning
tmp
align.sh
Definition: createJobs.py:716
double phi() const final
momentum azimuthal angle
void setP4(const LorentzVector &p4) final
set 4-momentum
edm::EDGetTokenT< reco::CaloMETCollection > caloMetBEToken_