ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/joshmt/PatTools/EmbedTree/src/EmbedTree.cc
Revision: 1.4
Committed: Thu Dec 6 21:59:26 2012 UTC (12 years, 4 months ago) by joshmt
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +85 -9 lines
Log Message:
used for v3 trees

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: EmbedTree
4 // Class: EmbedTree
5 //
6 /**\class EmbedTree EmbedTree.cc PatTools/EmbedTree/src/EmbedTree.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Joshua Thompson
15 // Created: Tue Nov 6 10:04:04 CST 2012
16 // $Id: EmbedTree.cc,v 1.2 2012/11/20 14:49:14 joshmt Exp $
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23
24 // user include files
25 #include "FWCore/Framework/interface/Frameworkfwd.h"
26 #include "FWCore/Framework/interface/EDAnalyzer.h"
27
28 #include "FWCore/Framework/interface/Event.h"
29 #include "FWCore/Framework/interface/MakerMacros.h"
30
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32
33 //my includes
34 #include "DataFormats/Common/interface/Handle.h"
35 #include "DataFormats/Common/interface/View.h"
36 #include "FWCore/ServiceRegistry/interface/Service.h"
37 #include "CommonTools/UtilAlgos/interface/TFileService.h"
38 #include "DataFormats/JetReco/interface/Jet.h"
39 #include "DataFormats/METReco/interface/MET.h"
40
41 #include "DataFormats/MuonReco/interface/Muon.h"
42 #include "DataFormats/PatCandidates/interface/Muon.h"
43 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
44
45 #include "DataFormats/VertexReco/interface/Vertex.h"
46
47 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
48
49 #include "DataFormats/Math/interface/deltaPhi.h"
50
51 //ROOT includes
52 #include "TTree.h"
53
54 //
55 // class declaration
56 //
57
58 class EmbedTree : public edm::EDAnalyzer {
59 public:
60 explicit EmbedTree(const edm::ParameterSet&);
61 ~EmbedTree();
62
63 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
64
65
66 private:
67 virtual void beginJob() ;
68 virtual void analyze(const edm::Event&, const edm::EventSetup&);
69 virtual void endJob() ;
70
71 virtual void beginRun(edm::Run const&, edm::EventSetup const&);
72 virtual void endRun(edm::Run const&, edm::EventSetup const&);
73 virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
74 virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
75
76 void fillTtbarDecayCode(const edm::Event&, const edm::EventSetup&);
77 void fillIsoTracks(const edm::Event&, const edm::EventSetup&);
78 void fillLeptons(const edm::Event&, const edm::EventSetup&);
79 void fillPreembedInfo(const edm::Event&, const edm::EventSetup&);
80
81 void clearTreeVariables();
82 int translateTopCode2PdgId(const int code);
83
84
85 // ----------member data ---------------------------
86
87 // -- config info
88 edm::InputTag vtxSrc_;
89 edm::InputTag jetSrc_, metSrc_;
90 edm::InputTag ttbarDecaySrc_;
91 edm::InputTag muonSrc_,electronSrc_;
92 edm::InputTag mindrSrc_;
93 edm::InputTag originalMuonSrc_;
94
95 // ntuple stuff
96
97 TTree* tree_;
98
99 int nGoodPV;
100
101 //jet counters
102 int njets70;
103 int njets50;
104 int njets30;
105
106 //TODO b tag info [can't until that is working in embedding]
107 //TODO top tagging, HepTopTagger
108 //TODO MT2
109 //TODO mT(b,MET)
110
111 //for embedded events only, quantities for the original muon
112 float minDRmuonJet;
113 int nseedmuons;
114 float seedmuonpt,seedmuoneta,seedmuonphi;
115 float seedWpt;
116
117 //MET, jet info
118 float DeltaPhiMetJet1;
119 float DeltaPhiMetJet2;
120 float DeltaPhiMetJet3;
121
122 //jet info
123 float jetpt1,jetphi1,jeteta1;
124
125 // float DeltaPhiMetIsotrack1; //can't do this until I have phi of isotrack
126 float DeltaPhiMetMuon1,DeltaPhiMetElectron1;
127
128 //MET
129 float MET;
130 float METphi;
131
132 //isolated tracks
133 float isotrackpt1;//,isotrackphi1, isotracketa1; (no dice here...would need to modify producer)
134
135 //leptons
136 float electronpt1,electroneta1,electronphi1;
137 float muonpt1,muoneta1,muonphi1;
138
139
140 //mc truth, for unembedded MC only
141 int ttbarDecayCode;
142 float genLeptonPt1,genLeptonPt2;
143 float genLeptonVisPt1,genLeptonVisPt2;
144 int genLeptonNProng1,genLeptonNProng2;
145 float genLeptonEta1,genLeptonEta2;
146 float genLeptonPhi1,genLeptonPhi2;
147 int genLeptonFlavor1,genLeptonFlavor2; //provides more info that the (legacy) ttbarDecayCode
148
149 };
150
151 void EmbedTree::clearTreeVariables() {
152
153 nGoodPV=0;
154
155 //to be incremented -- must be zero
156 njets70=0;
157 njets50=0;
158 njets30=0;
159
160 jetpt1=-1; jetphi1=-99; jeteta1=-99;
161
162 //generic bogus values
163 ttbarDecayCode=-1;
164 genLeptonPt1=-1; genLeptonPt2=-1;
165 genLeptonVisPt1=-1; genLeptonVisPt2=-1;
166 genLeptonNProng1=0; genLeptonNProng2=0;
167 genLeptonEta1=-99; genLeptonEta2=-99;
168 genLeptonPhi1=-99; genLeptonPhi2=-99;
169 genLeptonFlavor1=-1; genLeptonFlavor2=-1;
170
171 MET=-99;
172 METphi=-99;
173
174 isotrackpt1=-99; //isotrackphi1=-99; isotracketa1=-99;
175
176 electronpt1=-99; electroneta1=-99; electronphi1=-99;
177 muonpt1=-99; muoneta1=-99; muonphi1=-99;
178
179 DeltaPhiMetJet1=-99; DeltaPhiMetJet2=-99; DeltaPhiMetJet3=-99;
180 // DeltaPhiMetIsotrack1=-99;
181 DeltaPhiMetMuon1=-99;
182 DeltaPhiMetElectron1=-99;
183 minDRmuonJet=-99;
184 nseedmuons=0;
185 seedmuonpt=-99; seedmuoneta=-99; seedmuonphi=-99;
186 seedWpt=-99;
187
188 }
189
190 //
191 // constants, enums and typedefs
192 //
193
194 //
195 // static data member definitions
196 //
197
198 //
199 // constructors and destructor
200 //
201 EmbedTree::EmbedTree(const edm::ParameterSet& iConfig) :
202 vtxSrc_(iConfig.getParameter<edm::InputTag>("VtxSource")),
203 jetSrc_(iConfig.getParameter<edm::InputTag>("JetSource")),
204 metSrc_(iConfig.getParameter<edm::InputTag>("MetSource")),
205 ttbarDecaySrc_(iConfig.getParameter<edm::InputTag>("TtbarDecaySource")),
206 muonSrc_(iConfig.getParameter<edm::InputTag>("MuonSource")),
207 electronSrc_(iConfig.getParameter<edm::InputTag>("ElectronSource")),
208 mindrSrc_(iConfig.getParameter<edm::InputTag>("MinDRSource")),
209 originalMuonSrc_(iConfig.getParameter<edm::InputTag>("OriginalMuonSource")),
210 tree_(0)
211 {
212 //now do what ever initialization is needed
213
214 //now do what ever initialization is needed
215 edm::Service<TFileService> fs;
216 //branch creation for tree done in beginJob instead of here
217 tree_= fs->make<TTree>("reducedTree","reducedTree");
218
219
220 }
221
222
223 EmbedTree::~EmbedTree()
224 {
225
226 // do anything here that needs to be done at desctruction time
227 // (e.g. close files, deallocate resources etc.)
228
229 }
230
231
232 //
233 // member functions
234 //
235
236 // ------------ method called for each event ------------
237 void
238 EmbedTree::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
239 {
240 using namespace edm;
241
242 clearTreeVariables();
243
244 //important that MET is filled first
245 edm::Handle<edm::View<reco::MET> > mets;
246 iEvent.getByLabel(metSrc_, mets);
247 MET = mets->front().pt();
248 METphi = mets->front().phi();
249
250 //fill pv info
251 edm::Handle<edm::View<reco::Vertex> > pvs;
252 iEvent.getByLabel(vtxSrc_, pvs);
253 nGoodPV = pvs->size();
254
255 //do not use pat::Jet
256 edm::Handle<edm::View<reco::Jet> > jets;
257 iEvent.getByLabel(jetSrc_,jets);
258
259 for (edm::View<reco::Jet>::const_iterator jet = jets->begin(); jet != jets->end(); ++jet) {
260
261 if ( fabs(jet->eta()) < 2.4 ) {
262 float pt = jet->pt();
263 if ( pt >= 30) {
264 ++njets30;
265 if (njets30 ==1) DeltaPhiMetJet1 = std::abs(reco::deltaPhi(jet->phi(),METphi));
266 else if (njets30==2) DeltaPhiMetJet2 = std::abs(reco::deltaPhi(jet->phi(),METphi));
267 else if (njets30==3) DeltaPhiMetJet3 = std::abs(reco::deltaPhi(jet->phi(),METphi));
268 }
269 if (pt>=50) ++njets50;
270 if (pt>=70) ++njets70;
271
272 if (pt>jetpt1) { //lead jet info
273 jetpt1 = pt;
274 jeteta1 = jet->eta();
275 jetphi1 = jet->phi();
276 }
277 }
278 }
279
280 fillTtbarDecayCode(iEvent, iSetup);
281
282 fillIsoTracks(iEvent,iSetup);
283
284 fillLeptons(iEvent,iSetup);
285
286 fillPreembedInfo(iEvent,iSetup);
287
288 //fill tree
289 tree_->Fill();
290 }
291
292 void EmbedTree::fillPreembedInfo(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
293
294 if ( std::string(mindrSrc_.label()).empty()) minDRmuonJet = -1;
295 else {
296 edm::Handle<float > mdr;
297 iEvent.getByLabel(mindrSrc_.label(),"minDeltaRToJet",mdr);
298 minDRmuonJet = *mdr;
299 }
300
301 if (std::string(originalMuonSrc_.label()).empty()) {
302 //don't really need to do anything. the variables are reset in clearTreeVariables()
303 }
304 else {
305 edm::Handle<edm::View<reco::PFCandidate> > mulist;
306 iEvent.getByLabel(originalMuonSrc_,mulist);
307
308 edm::Handle<edm::View<pat::Muon> > mulistpat;
309 if (mulist.failedToGet()) iEvent.getByLabel(originalMuonSrc_,mulistpat);
310
311 nseedmuons=0;
312 if (mulist.failedToGet() ) { //use pat muons
313 // std::cout<<" [using PAT muons]"<<std::endl;
314 for (edm::View<pat::Muon>::const_iterator themu = mulistpat->begin(); themu != mulistpat->end(); ++themu) {
315 ++nseedmuons;
316 if(nseedmuons==1) { //store info for the first guy only
317 seedmuonpt = themu->pt();
318 seedmuoneta = themu->eta();
319 seedmuonphi = themu->phi();
320 }
321 }
322 }
323 else {//use pf cands
324 // std::cout<<" [using pf muons]"<<std::endl;
325 for (edm::View<reco::PFCandidate>::const_iterator themu = mulist->begin(); themu != mulist->end(); ++themu) {
326 ++nseedmuons;
327 if(nseedmuons==1) { //store info for the first guy only
328 seedmuonpt = themu->pt();
329 seedmuoneta = themu->eta();
330 seedmuonphi = themu->phi();
331 }
332 }
333 }
334 }
335
336 //vector<reco::CompositeCandidate> "WtoMuNu" "" "EMBED"
337 edm::Handle<edm::View<reco::CompositeCandidate> > wtomunu;
338 iEvent.getByLabel("WtoMuNu",wtomunu);
339 if (!wtomunu.failedToGet()) {
340 int nn=0;
341 for (edm::View<reco::CompositeCandidate>::const_iterator thew = wtomunu->begin(); thew != wtomunu->end(); ++thew) {
342 if (nn==0) {
343 seedWpt = thew->pt();
344 nn++;
345 }
346 }
347 }
348
349 }
350
351 void EmbedTree::fillIsoTracks(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
352
353 //for now the only thing we store is the pT of the leading isolated track
354
355 edm::Handle<std::vector<float> > cand_dzpv;
356 iEvent.getByLabel("trackIsolationMaker","pfcandsdzpv",cand_dzpv);
357
358 edm::Handle<std::vector<float> > cand_pt;
359 iEvent.getByLabel("trackIsolationMaker","pfcandspt",cand_pt);
360
361 edm::Handle<std::vector<float> > cand_trkiso;
362 iEvent.getByLabel("trackIsolationMaker","pfcandstrkiso",cand_trkiso);
363
364 edm::Handle<std::vector<int> > cand_chg;
365 iEvent.getByLabel("trackIsolationMaker","pfcandschg",cand_chg);
366
367 size_t n = cand_dzpv->size();
368 for ( size_t ii = 0; ii<n; ii++) {
369
370 if ( cand_chg->at(ii) != 0 && (fabs(cand_dzpv->at(ii)) < 0.05)) { //charged and compatible with PV
371
372 //No pT cut! always store lead iso track
373 //if (cand_pt->at(ii) > 10) { //pT >10
374
375 if (cand_trkiso->at(ii)/cand_pt->at(ii) <0.10) {
376 if ( cand_pt->at(ii) > isotrackpt1) {
377 isotrackpt1 = cand_pt->at(ii);
378 }
379 }
380 // }
381 }
382 }
383
384 }
385
386 int EmbedTree::translateTopCode2PdgId(const int code) {
387
388
389 //undo my old scheme and return the well-known pdg id
390
391 if (code==2) return 11; //e
392 else if (code==3) return 13;//mu
393 else if (code>=4) return 15;//tau
394
395 return 0;
396
397 }
398
399 void EmbedTree::fillTtbarDecayCode(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
400
401 if ( std::string(ttbarDecaySrc_.label()).empty()) ttbarDecayCode = -1;
402 else {
403 edm::Handle<int > code;
404 iEvent.getByLabel(ttbarDecaySrc_.label(),"ttbarDecayCode",code);
405 ttbarDecayCode = *code;
406
407 //also get the generator-level lepton pT and eta
408 edm::Handle<std::vector<float> > genleppt;
409 edm::Handle<std::vector<float> > gentauvispt;
410 edm::Handle<std::vector<int> > gentaunprong;
411 edm::Handle<std::vector<float> > genlepeta;
412 edm::Handle<std::vector<float> > genlepphi;
413 edm::Handle<std::vector<int> > genlepflavor;
414 iEvent.getByLabel(ttbarDecaySrc_.label(),"lepGenPt",genleppt);
415 iEvent.getByLabel(ttbarDecaySrc_.label(),"tauGenVisPt",gentauvispt);
416 iEvent.getByLabel(ttbarDecaySrc_.label(),"tauGenNProng",gentaunprong);
417 iEvent.getByLabel(ttbarDecaySrc_.label(),"lepGenPhi",genlepphi);
418 iEvent.getByLabel(ttbarDecaySrc_.label(),"lepGenEta",genlepeta);
419 iEvent.getByLabel(ttbarDecaySrc_.label(),"topDecayCode",genlepflavor);
420
421 //sort by pT, so that the highest pT lepton will be in the *1 variables
422 int i=0; int j=1;
423 if ( (*genleppt)[0] > (*genleppt)[1] ) {
424 i=0;
425 j=1;
426 }
427 else {
428 i=1;
429 j=0;
430 }
431
432 genLeptonPt1 = (*genleppt)[i];
433 genLeptonPt2 = (*genleppt)[j];
434
435 genLeptonVisPt1 = (*gentauvispt)[i];
436 genLeptonVisPt2 = (*gentauvispt)[j];
437
438 genLeptonNProng1 = (*gentaunprong)[i];
439 genLeptonNProng2 = (*gentaunprong)[j];
440
441 genLeptonEta1 = (*genlepeta)[i];
442 genLeptonEta2 = (*genlepeta)[j];
443
444 genLeptonPhi1 = (*genlepphi)[i];
445 genLeptonPhi2 = (*genlepphi)[j];
446
447 genLeptonFlavor1 = translateTopCode2PdgId((*genlepflavor)[i]);
448 genLeptonFlavor2 = translateTopCode2PdgId((*genlepflavor)[j]);
449 }
450
451 }
452
453
454 void EmbedTree::fillLeptons(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
455
456 //fill info on highest pt leptons
457
458 //i'll bet a templated function would serve well here, but let's keep it simple for now
459
460 edm::Handle<edm::View<reco::GsfElectron> > electrons;
461 iEvent.getByLabel(electronSrc_, electrons);
462 for (edm::View<reco::GsfElectron>::const_iterator ie = electrons->begin(); ie != electrons->end(); ++ie) {
463 if ( ie->et() > electronpt1) {
464 electronpt1 = ie->et();
465 electronphi1 = ie->phi();
466 electroneta1 = ie->eta();
467 DeltaPhiMetElectron1 = std::abs(reco::deltaPhi(electronphi1,METphi));
468 }
469 }
470
471
472 edm::Handle<edm::View<reco::Muon> > muons;
473 iEvent.getByLabel(muonSrc_, muons);
474 for (edm::View<reco::Muon >::const_iterator im = muons->begin(); im != muons->end(); ++im) {
475 if ( im->pt() > muonpt1) {
476 muonpt1 = im->pt();
477 muonphi1 = im->phi();
478 muoneta1 = im->eta();
479 DeltaPhiMetMuon1 = std::abs(reco::deltaPhi(muonphi1,METphi));
480 }
481 }
482
483
484 }
485
486
487 // ------------ method called once each job just before starting event loop ------------
488
489 void
490 EmbedTree::beginJob()
491 {
492
493 tree_->Branch("nGoodPV",&nGoodPV,"nGoodPV/I");
494
495 tree_->Branch("njets70",&njets70,"njets70/I");
496 tree_->Branch("njets50",&njets50,"njets50/I");
497 tree_->Branch("njets30",&njets30,"njets30/I");
498
499 tree_->Branch("MET",&MET,"MET/F");
500 tree_->Branch("METphi",&METphi,"METphi/F");
501
502 tree_->Branch("DeltaPhiMetJet1",&DeltaPhiMetJet1,"DeltaPhiMetJet1/F");
503 tree_->Branch("DeltaPhiMetJet2",&DeltaPhiMetJet2,"DeltaPhiMetJet2/F");
504 tree_->Branch("DeltaPhiMetJet3",&DeltaPhiMetJet3,"DeltaPhiMetJet3/F");
505
506 tree_->Branch("DeltaPhiMetMuon1",&DeltaPhiMetMuon1,"DeltaPhiMetMuon1/F");
507 tree_->Branch("DeltaPhiMetElectron1",&DeltaPhiMetElectron1,"DeltaPhiMetElectron1/F");
508
509 tree_->Branch("jetpt1",&jetpt1,"jetpt1/F");
510 tree_->Branch("jetphi1",&jetphi1,"jetphi1/F");
511 tree_->Branch("jeteta1",&jeteta1,"jeteta1/F");
512
513 tree_->Branch("isotrackpt1",&isotrackpt1,"isotrackpt1/F");
514
515 tree_->Branch("electronpt1",&electronpt1,"electronpt1/F");
516 tree_->Branch("electronphi1",&electronphi1,"electronphi1/F");
517 tree_->Branch("electroneta1",&electroneta1,"electroneta1/F");
518
519 tree_->Branch("muonpt1",&muonpt1,"muonpt1/F");
520 tree_->Branch("muonphi1",&muonphi1,"muonphi1/F");
521 tree_->Branch("muoneta1",&muoneta1,"muoneta1/F");
522
523 tree_->Branch("nseedmuons",&nseedmuons,"nseedmuons/I");
524 tree_->Branch("seedmuonpt",&seedmuonpt,"seedmuonpt/F");
525 tree_->Branch("seedmuonphi",&seedmuonphi,"seedmuonphi/F");
526 tree_->Branch("seedmuoneta",&seedmuoneta,"seedmuoneta/F");
527 tree_->Branch("seedWpt",&seedWpt,"seedWpt/F");
528
529 tree_->Branch("minDRmuonJet",&minDRmuonJet,"minDRmuonJet/F");
530
531
532 tree_->Branch("ttbarDecayCode",&ttbarDecayCode,"ttbarDecayCode/I");
533 tree_->Branch("genLeptonFlavor1",&genLeptonFlavor1,"genLeptonFlavor1/I");
534 tree_->Branch("genLeptonPt1",&genLeptonPt1,"genLeptonPt1/F");
535 tree_->Branch("genLeptonVisPt1",&genLeptonVisPt1,"genLeptonVisPt1/F");
536 tree_->Branch("genLeptonNProng1",&genLeptonNProng1,"genLeptonNProng1/I");
537 tree_->Branch("genLeptonEta1",&genLeptonEta1,"genLeptonEta1/F");
538 tree_->Branch("genLeptonPhi1",&genLeptonPhi1,"genLeptonPhi1/F");
539
540 tree_->Branch("genLeptonFlavor2",&genLeptonFlavor2,"genLeptonFlavor2/I");
541 tree_->Branch("genLeptonPt2",&genLeptonPt2,"genLeptonPt2/F");
542 tree_->Branch("genLeptonVisPt2",&genLeptonVisPt2,"genLeptonVisPt2/F");
543 tree_->Branch("genLeptonNProng2",&genLeptonNProng2,"genLeptonNProng2/I");
544 tree_->Branch("genLeptonEta2",&genLeptonEta2,"genLeptonEta2/F");
545 tree_->Branch("genLeptonPhi2",&genLeptonPhi2,"genLeptonPhi2/F");
546
547
548 }
549
550 // ------------ method called once each job just after ending the event loop ------------
551 void
552 EmbedTree::endJob()
553 {
554 }
555
556 // ------------ method called when starting to processes a run ------------
557 void
558 EmbedTree::beginRun(edm::Run const&, edm::EventSetup const&)
559 {
560 }
561
562 // ------------ method called when ending the processing of a run ------------
563 void
564 EmbedTree::endRun(edm::Run const&, edm::EventSetup const&)
565 {
566 }
567
568 // ------------ method called when starting to processes a luminosity block ------------
569 void
570 EmbedTree::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
571 {
572 }
573
574 // ------------ method called when ending the processing of a luminosity block ------------
575 void
576 EmbedTree::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
577 {
578 }
579
580 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
581 void
582 EmbedTree::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
583 //The following says we do not know what parameters are allowed so do no validation
584 // Please change this to state exactly what you do use, even if it is no parameters
585 edm::ParameterSetDescription desc;
586 desc.setUnknown();
587 descriptions.addDefault(desc);
588 }
589
590 //define this as a plug-in
591 DEFINE_FWK_MODULE(EmbedTree);