test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
L1JetRecoTreeProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: L1Trigger/L1TNtuples
4 // Class: L1JetRecoTreeProducer
5 //
14 // system include files
15 #include <memory>
16 
17 // framework
26 
27 // cond formats
31 
32 // data formats
40 
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 
59 public:
62 
63 
64 private:
65  virtual void beginJob(void) ;
66  virtual void analyze(const edm::Event&, const edm::EventSetup&);
67  virtual void endJob();
68 
74 
76 
77  bool jetId(const reco::PFJet& jet);
78 
79 public:
82 
83 private:
84 
85  // output file
87 
88  // tree
89  TTree * tree_;
90 
91  // EDM input tags
93  // edm::EDGetTokenT<reco::CaloJetCollection> caloJetToken_;
94  // edm::EDGetTokenT<edm::ValueMap<reco::JetID> > caloJetIdToken_;
96 
100 
101 
102  // debug stuff
105  double jetetaMax_;
106  unsigned int maxCl_;
107  unsigned int maxJet_;
108  unsigned int maxVtx_;
109  unsigned int maxTrk_;
110 
115 
116 };
117 
118 
120  pfJetsMissing_(false),
121  pfMetMissing_(false),
122  pfJetCorrMissing_(false),
123  caloMetMissing_(false),
124  caloMetBEMissing_(false)
125 {
126 
127  // caloJetToken_ = consumes<reco::CaloJetCollection>(iConfig.getUntrackedParameter("caloJetToken",edm::InputTag("ak4CaloJets")));
128  pfJetToken_ = consumes<reco::PFJetCollection>(iConfig.getUntrackedParameter("pfJetToken",edm::InputTag("ak4PFJetsCHS")));
129  // caloJetIdToken_ = consumes<edm::ValueMap<reco::JetID> >(iConfig.getUntrackedParameter("jetIdToken",edm::InputTag("ak4JetID")));
130  jecToken_ = consumes<reco::JetCorrector>(iConfig.getUntrackedParameter<edm::InputTag>("jecToken",edm::InputTag("ak4PFCHSL1FastL2L3ResidualCorrector")));
131 
132  pfMetToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("pfMetToken",edm::InputTag("pfMet")));
133  caloMetToken_ = consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetToken",edm::InputTag("caloMet")));
134  caloMetBEToken_ = consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetBEToken",edm::InputTag("caloMetBE")));
135 
136  jetptThreshold_ = iConfig.getParameter<double> ("jetptThreshold");
137  jetetaMax_ = iConfig.getParameter<double> ("jetetaMax");
138  maxJet_ = iConfig.getParameter<unsigned int>("maxJet");
139 
142 
143  // set up output
144  tree_=fs_->make<TTree>("JetRecoTree", "JetRecoTree");
145  tree_->Branch("Jet", "L1Analysis::L1AnalysisRecoJetDataFormat", &jet_data, 32000, 3);
146  tree_->Branch("Sums", "L1Analysis::L1AnalysisRecoMetDataFormat", &met_data, 32000, 3);
147 
148 }
149 
150 
152 {
153 
154  // do anything here that needs to be done at desctruction time
155  // (e.g. close files, deallocate resources etc.)
156 
157 }
158 
159 
160 //
161 // member functions
162 //
163 
164 // ------------ method called to for each event ------------
166 {
167 
168  jet_data->Reset();
169  met_data->Reset();
170 
171  // get jets
173  iEvent.getByToken(pfJetToken_, pfJets);
174 
175  //get sums
177  iEvent.getByToken(pfMetToken_, pfMet);
178 
179  // get jet ID
180  // edm::Handle<edm::ValueMap<reco::JetID> > jetsID;
181  //iEvent.getByLabel(jetIdTag_,jetsID);
182 
184  iEvent.getByToken(jecToken_, pfJetCorr);
185 
187  iEvent.getByToken(caloMetToken_, caloMet);
188 
190  iEvent.getByToken(caloMetBEToken_, caloMetBE);
191 
192  if (pfJets.isValid()) {
193 
194  jet_data->nJets=0;
195 
196  doPFJets(pfJets);
197 
198  }
199  else {
200  if (!pfJetsMissing_) {edm::LogWarning("MissingProduct") << "PFJets not found. Branch will not be filled" << std::endl;}
201  pfJetsMissing_ = true;
202  }
203 
204  if (pfJetCorr.isValid()) {
205 
206  doPFJetCorr(pfJets,pfJetCorr);
207 
208  }
209  else {
210  if (!pfJetCorrMissing_) {edm::LogWarning("MissingProduct") << "Jet Corrector not found. Branch will not be filled" << std::endl;}
211  pfJetCorrMissing_ = true;
212  }
213 
214  if (pfMet.isValid()) {
215 
216  doPFMet(pfMet);
217 
218  }
219  else {
220  if (!pfMetMissing_) {edm::LogWarning("MissingProduct") << "PFMet not found. Branch will not be filled" << std::endl;}
221  pfMetMissing_ = true;
222  }
223 
224  if (caloMet.isValid()) {
225 
226  doCaloMet(caloMet);
227 
228  }
229  else {
230  if (!caloMetMissing_) {edm::LogWarning("MissingProduct") << "CaloMet not found. Branch will not be filled" << std::endl;}
231  caloMetMissing_ = true;
232  }
233 
234  if (caloMetBE.isValid()) {
235 
236  doCaloMetBE(caloMetBE);
237 
238  }
239  else {
240  if (!caloMetBEMissing_) {edm::LogWarning("MissingProduct") << "CaloMetBE not found. Branch will not be filled" << std::endl;}
241  caloMetBEMissing_ = true;
242  }
243 
244  tree_->Fill();
245 
246 }
247 
248 
249 void
251 
252 
253  for( auto it=caloJets->begin();
254  it!=caloJets->end() && jet_data->nJets < maxJet_;
255  ++it) {
256 
257  jet_data->et.push_back(it->et());
258  jet_data->eta.push_back(it->eta());
259  jet_data->phi.push_back(it->phi());
260  jet_data->e.push_back(it->energy());
261  jet_data->isPF.push_back(false);
262 
263  jet_data->eEMF.push_back(it->emEnergyFraction());
264  jet_data->eEmEB.push_back(it->emEnergyInEB());
265  jet_data->eEmEE.push_back(it->emEnergyInEE());
266  jet_data->eEmHF.push_back(it->emEnergyInHF());
267  jet_data->eHadHB.push_back(it->hadEnergyInHB());
268  jet_data->eHadHE.push_back(it->hadEnergyInHE());
269  jet_data->eHadHO.push_back(it->hadEnergyInHO());
270  jet_data->eHadHF.push_back(it->hadEnergyInHF());
271  jet_data->eMaxEcalTow.push_back(it->maxEInEmTowers());
272  jet_data->eMaxHcalTow.push_back(it->maxEInHadTowers());
273  jet_data->towerArea.push_back(it->towersArea());
274 
275  jet_data->nJets++;
276 
277  }
278 
279 
280 }
281 
282 
283 void
285 
286 
287  for( auto it=pfJets->begin();
288  it!=pfJets->end() && jet_data->nJets < maxJet_;
289  ++it) {
290  jet_data->et.push_back(it->et());
291  jet_data->eta.push_back(it->eta());
292  jet_data->phi.push_back(it->phi());
293  jet_data->e.push_back(it->energy());
294  jet_data->isPF.push_back(true);
295 
296  jet_data->chef.push_back(it->chargedHadronEnergyFraction());
297  jet_data->nhef.push_back(it->neutralHadronEnergyFraction());
298  jet_data->pef.push_back(it->photonEnergyFraction());
299  jet_data->eef.push_back(it->electronEnergyFraction());
300  jet_data->mef.push_back(it->muonEnergyFraction());
301  jet_data->hfhef.push_back(it->HFHadronEnergyFraction());
302  jet_data->hfemef.push_back(it->HFEMEnergyFraction());
303  jet_data->chMult.push_back(it->chargedHadronMultiplicity());
304  jet_data->nhMult.push_back(it->neutralHadronMultiplicity());
305  jet_data->phMult.push_back(it->photonMultiplicity());
306  jet_data->elMult.push_back(it->electronMultiplicity());
307  jet_data->muMult.push_back(it->muonMultiplicity());
308  jet_data->hfhMult.push_back(it->HFHadronMultiplicity());
309  jet_data->hfemMult.push_back(it->HFEMMultiplicity());
310 
311  jet_data->cemef.push_back(it->chargedEmEnergyFraction());
312  jet_data->cmef.push_back(it->chargedMuEnergyFraction());
313  jet_data->nemef.push_back(it->neutralEmEnergyFraction());
314  jet_data->cMult.push_back(it->chargedMultiplicity());
315  jet_data->nMult.push_back(it->neutralMultiplicity());
316 
317  jet_data->nJets++;
318 
319  }
320 
321 }
322 
323 
324 void
326 
327 
328  float corrFactor = 1.;
329  uint nJets = 0;
330 
331  float mHx = 0;
332  float mHy = 0;
333 
334  met_data->Ht = 0;
335  met_data->mHt = -999.;
336  met_data->mHtPhi = -999.;
337 
338 
339  for( auto it=pfJets->begin();
340  it!=pfJets->end() && nJets < maxJet_;
341  ++it) {
342 
343  corrFactor = pfJetCorr.product()->correction(*it);
344 
345  jet_data->etCorr.push_back(it->et()*corrFactor);
346  jet_data->corrFactor.push_back(corrFactor);
347 
348  nJets++;
349 
350  if (it->pt()*corrFactor > jetptThreshold_ && fabs(it->eta())<jetetaMax_) {
351  mHx += -1.*it->px()*corrFactor;
352  mHy += -1.*it->py()*corrFactor;
353  met_data->Ht += it->pt()*corrFactor;
354  }
355 
356  }
357 
358  TVector2 *tv2 = new TVector2(mHx,mHy);
359  met_data->mHt = tv2->Mod();
360  met_data->mHtPhi = tv2->Phi();
361 
362  // std::vector< std::pair<float,float> > corrJetEtsAndCorrs;
363 
364  // //get jet correction and fill corrected jet ets and corrections
365  // for( auto it=pfJets->begin(); it!=pfJets->end(); ++it)
366  // {
367  // float corr = pfJetCorr.product()->correction(*it);
368  // std::pair<float,float> corrJetEtAndCorr(corr*it->et(),corr);
369  // corrJetEtsAndCorrs.push_back(corrJetEtAndCorr);
370  // }
371 
372  // // sort corrected jet ets and correction factors
373  // // by corrected jet et
374  // std::sort(corrJetEtsAndCorrs.rbegin(),corrJetEtsAndCorrs.rend());
375 
376  // //fill jet data array with sorted jet ets and corr factors
377  // std::vector<std::pair<float,float> >::iterator it;
378  // uint nJets = 0;
379 
380  // for(it = corrJetEtsAndCorrs.begin(); it != corrJetEtsAndCorrs.end() && nJets < maxJet_; ++it){
381  // jet_data->etCorr.push_back(it->first);
382  // jet_data->corrFactor.push_back(it->second);
383  // nJets++
384  // }
385 
386 }
387 
388 void
390 
391  const reco::PFMETCollection *metCol = pfMet.product();
392  const reco::PFMET theMet = metCol->front();
393 
394  met_data->met = theMet.et();
395  met_data->metPhi = theMet.phi();
396  met_data->sumEt = theMet.sumEt();
397 
398 }
399 
400 void
402 
403  const reco::CaloMETCollection *metCol = caloMet.product();
404  const reco::CaloMET theMet = metCol->front();
405 
406  met_data->caloMet = theMet.et();
407  met_data->caloMetPhi = theMet.phi();
408  met_data->caloSumEt = theMet.sumEt();
409 
410 }
411 
412 void
414 
415  const reco::CaloMETCollection *metCol = caloMetBE.product();
416  const reco::CaloMET theMet = metCol->front();
417 
418  met_data->caloMetBE = theMet.et();
419  met_data->caloMetPhiBE = theMet.phi();
420  met_data->caloSumEtBE = theMet.sumEt();
421 
422 }
423 
424 bool
426 
427  bool tmp = true;
428 
429  tmp &= jet.neutralHadronEnergyFraction() < 0.9 ;
430  tmp &= jet.neutralEmEnergyFraction() < 0.9 ;
431  tmp &= (jet.chargedMultiplicity() + jet.neutralMultiplicity()) > 1 ;
432  tmp &= jet.muonEnergyFraction() < 0.8 ;
433  if (fabs(jet.eta()) < 2.4) {
434  tmp &= jet.chargedHadronEnergyFraction() > 0.0 ;
435  tmp &= jet.chargedMultiplicity() > 0 ;
436  tmp &= jet.chargedEmEnergyFraction() < 0.9 ;
437  }
438  if (fabs(jet.eta()) > 3.0) {
439  tmp &= jet.neutralEmEnergyFraction() < 0.9 ;
440  tmp &= jet.neutralMultiplicity() > 10 ;
441  }
442 
443  // our custom selection
444  tmp &= jet.muonMultiplicity() == 0;
445  tmp &= jet.electronMultiplicity() == 0;
446 
447  return tmp;
448 
449 }
450 
451 
452 // ------------ method called once each job just before starting event loop ------------
453 void
455 {
456 }
457 
458 // ------------ method called once each job just after ending the event loop ------------
459 void
461 }
462 
463 //define this as a plug-in
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
void doPFMet(edm::Handle< reco::PFMETCollection > pfMet)
L1Analysis::L1AnalysisRecoMetDataFormat * met_data
tuple caloMet
____________________________________________________________________________||
Definition: CaloMET_cfi.py:4
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
edm::Service< TFileService > fs_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
float chargedHadronEnergyFraction() const
chargedHadronEnergyFraction
Definition: PFJet.h:100
void doPFJetCorr(edm::Handle< reco::PFJetCollection > pfJets, edm::Handle< reco::JetCorrector > pfJetCorr)
edm::EDGetTokenT< reco::CaloMETCollection > caloMetToken_
virtual double phi() const final
momentum azimuthal angle
void doPFJets(edm::Handle< reco::PFJetCollection > pfJets)
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
void doCaloMetBE(edm::Handle< reco::CaloMETCollection > caloMetBE)
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:155
Jets made from PFObjects.
Definition: PFJet.h:21
L1Analysis::L1AnalysisRecoJetDataFormat * jet_data
float neutralEmEnergyFraction() const
neutralEmEnergyFraction
Definition: PFJet.h:152
int iEvent
Definition: GenABIO.cc:230
double sumEt() const
Definition: MET.h:56
edm::EDGetTokenT< reco::JetCorrector > jecToken_
tuple pfMet
____________________________________________________________________________||
Definition: PFMET_cfi.py:4
float neutralHadronEnergyFraction() const
neutralHadronEnergyFraction
Definition: PFJet.h:104
edm::EDGetTokenT< reco::PFJetCollection > pfJetToken_
bool isValid() const
Definition: HandleBase.h:75
virtual void analyze(const edm::Event &, const edm::EventSetup &)
void doCaloMet(edm::Handle< reco::CaloMETCollection > caloMet)
int neutralMultiplicity() const
neutralMultiplicity
Definition: PFJet.h:157
edm::EDGetTokenT< reco::PFMETCollection > pfMetToken_
float chargedEmEnergyFraction() const
chargedEmEnergyFraction
Definition: PFJet.h:144
std::vector< reco::CaloMET > CaloMETCollection
collection of CaloMET objects
T const * product() const
Definition: Handle.h:81
L1JetRecoTreeProducer(const edm::ParameterSet &)
void doCaloJets(edm::Handle< reco::CaloJetCollection > caloJets)
int muonMultiplicity() const
muonMultiplicity
Definition: PFJet.h:135
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
std::vector< reco::PFMET > PFMETCollection
collection of PFMET objects
bool jetId(const reco::PFJet &jet)
virtual double et() const final
transverse energy
virtual double eta() const final
momentum pseudorapidity
volatile std::atomic< bool > shutdown_flag false
float muonEnergyFraction() const
muonEnergyFraction
Definition: PFJet.h:116
tuple pfJets
Definition: pfJets_cff.py:8
int electronMultiplicity() const
electronMultiplicity
Definition: PFJet.h:133
edm::EDGetTokenT< reco::CaloMETCollection > caloMetBEToken_
tuple caloMetBE
____________________________________________________________________________||
Definition: CaloMET_cfi.py:19