ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
Revision: 1.13
Committed: Thu Apr 19 14:47:13 2012 UTC (13 years ago) by peiffer
Content type: text/plain
Branch: MAIN
Changes since 1.12: +3 -1 lines
Log Message:
isolation

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: NtupleWriter
4 // Class: NtupleWriter
5 //
6 /**\class NtupleWriter NtupleWriter.cc NtupleWriter/src/NtupleWriter.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Thomas Peiffer,,,Uni Hamburg
15 // Created: Tue Mar 13 08:43:34 CET 2012
16 // $Id: NtupleWriter.cc,v 1.12 2012/04/19 12:03:28 peiffer Exp $
17 //
18 //
19
20 #include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21
22
23
24 //
25 // constants, enums and typedefs
26 //
27
28 //
29 // static data member definitions
30 //
31
32 //
33 // constructors and destructor
34 //
35 NtupleWriter::NtupleWriter(const edm::ParameterSet& iConfig)
36
37 {
38 //now do what ever initialization is needed
39 //edm::Service<TFileService> fs;
40 //tr = fs->make<TTree>("AnalysisTree","AnalysisTree");
41
42 fileName = iConfig.getParameter<std::string>("fileName");
43 outfile = new TFile(fileName,"RECREATE");
44 outfile->cd();
45 tr = new TTree("AnalysisTree","AnalysisTree");
46
47 std::string name;
48
49 doElectrons = iConfig.getParameter<bool>("doElectrons");
50 doMuons = iConfig.getParameter<bool>("doMuons");
51 doTaus = iConfig.getParameter<bool>("doTaus");
52 doJets = iConfig.getParameter<bool>("doJets");
53 doPhotons = iConfig.getParameter<bool>("doPhotons");
54 doMET = iConfig.getParameter<bool>("doMET");
55 doGenInfo = iConfig.getParameter<bool>("doGenInfo");
56 doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
57 doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
58 doPV = iConfig.getParameter<bool>("doPV");
59 doTopJets = iConfig.getParameter<bool>("doTopJets");
60 doTrigger = iConfig.getParameter<bool>("doTrigger");
61
62 // initialization of tree variables
63
64 tr->Branch("run",&run);
65 tr->Branch("event",&event);
66 tr->Branch("luminosityBlock",&luminosityBlock);
67 tr->Branch("isRealData",&isRealData);
68 tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
69 if(doLumiInfo){
70 tr->Branch("intgRecLumi",&intgRecLumi);
71 tr->Branch("intgDelLumi",&intgDelLumi);
72 }
73 if(doPV){
74 tr->Branch("beamspot_x0",&beamspot_x0);
75 tr->Branch("beamspot_y0",&beamspot_y0);
76 tr->Branch("beamspot_z0",&beamspot_z0);
77 }
78 if(doElectrons){
79 electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
80 for(size_t j=0; j< electron_sources.size(); ++j){
81 tr->Branch( electron_sources[j].c_str(), "std::vector<Electron>", &eles[j]);
82 }
83 }
84 if(doMuons){
85 muon_sources = iConfig.getParameter<std::vector<std::string> >("muon_sources");
86 for(size_t j=0; j< muon_sources.size(); ++j){
87 tr->Branch( muon_sources[j].c_str(), "std::vector<Muon>", &mus[j]);
88 }
89 }
90 if(doTaus){
91 tau_sources = iConfig.getParameter<std::vector<std::string> >("tau_sources");
92 tau_ptmin = iConfig.getParameter<double> ("tau_ptmin");
93 tau_etamax = iConfig.getParameter<double> ("tau_etamax");
94 for(size_t j=0; j< tau_sources.size(); ++j){
95 tr->Branch( tau_sources[j].c_str(), "std::vector<Tau>", &taus[j]);
96 }
97 }
98 if(doJets){
99 jet_sources = iConfig.getParameter<std::vector<std::string> >("jet_sources");
100 jet_ptmin = iConfig.getParameter<double> ("jet_ptmin");
101 jet_etamax = iConfig.getParameter<double> ("jet_etamax");
102 for(size_t j=0; j< jet_sources.size(); ++j){
103 tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
104 }
105 }
106 if(doTopJets){
107 topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
108 topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
109 topjet_etamax = iConfig.getParameter<double> ("topjet_etamax");
110 for(size_t j=0; j< topjet_sources.size(); ++j){
111 tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
112 }
113 }
114 if(doPhotons){
115 photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
116 for(size_t j=0; j< photon_sources.size(); ++j){
117 tr->Branch( photon_sources[j].c_str(), "std::vector<Photon>", &phs[j]);
118 }
119 }
120 if(doMET){
121 met_sources = iConfig.getParameter<std::vector<std::string> >("met_sources");
122 for(size_t j=0; j< met_sources.size(); ++j){
123 tr->Branch(met_sources[j].c_str(), "MET", &met[j]);
124 }
125 }
126 if(doPV){
127 pv_sources = iConfig.getParameter<std::vector<std::string> >("pv_sources");
128 for(size_t j=0; j< pv_sources.size(); ++j){
129 tr->Branch( pv_sources[j].c_str(), "std::vector<PrimaryVertex>", &pvs[j]);
130 }
131 }
132 if(doGenInfo){
133 tr->Branch("genInfo","GenInfo",&genInfo);
134 tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
135 }
136 if(doTrigger){
137 trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
138 //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
139 tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
140 tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
141 tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
142 tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
143 }
144 newrun = true;
145 }
146
147
148 NtupleWriter::~NtupleWriter()
149 {
150
151 // do anything here that needs to be done at desctruction time
152 // (e.g. close files, deallocate resources etc.)
153
154 }
155
156
157 //
158 // member functions
159 //
160
161 // ------------ method called for each event ------------
162 void
163 NtupleWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
164 {
165 // using namespace edm;
166
167
168 // ------------- common variables -------------
169
170 run = iEvent.id().run();
171 event = iEvent.id().event();
172 luminosityBlock = iEvent.luminosityBlock();
173 isRealData = iEvent.isRealData();
174
175 if(isRealData){
176 edm::Handle<bool> bool_handle;
177 iEvent.getByLabel(edm::InputTag("HBHENoiseFilterResultProducer","HBHENoiseFilterResult"),bool_handle);
178 HBHENoiseFilterResult = *(bool_handle.product());
179 }
180 else HBHENoiseFilterResult = false;
181
182 // ------------- primary vertices and beamspot -------------
183
184 if(doPV){
185 for(size_t j=0; j< pv_sources.size(); ++j){
186 pvs[j].clear();
187
188 edm::Handle< std::vector<reco::Vertex> > pv_handle;
189 iEvent.getByLabel(pv_sources[j], pv_handle);
190 const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
191
192 for (unsigned int i = 0; i < reco_pvs.size(); ++i) {
193 reco::Vertex reco_pv = reco_pvs[i];
194
195 PrimaryVertex pv;
196 pv.x = reco_pv.x();
197 pv.y = reco_pv.y();
198 pv.z = reco_pv.z();
199 pv.nTracks = reco_pv.nTracks();
200 //pv.isValid = reco_pv.isValid();
201 pv.chi2 = reco_pv.chi2();
202 pv.ndof = reco_pv.ndof();
203
204 pvs[j].push_back(pv);
205 }
206 }
207
208 edm::Handle<reco::BeamSpot> beamSpot;
209 iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
210 const reco::BeamSpot & bsp = *beamSpot;
211
212 beamspot_x0 = bsp.x0();
213 beamspot_y0 = bsp.y0();
214 beamspot_z0 = bsp.z0();
215 }
216
217 // ------------- electrons -------------
218 if(doElectrons){
219
220 edm::Handle<reco::ConversionCollection> hConversions;
221 iEvent.getByLabel("allConversions", hConversions);
222
223 edm::Handle<reco::BeamSpot> beamSpot;
224 iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
225 const reco::BeamSpot & bsp = *beamSpot;
226
227 for(size_t j=0; j< electron_sources.size(); ++j){
228 eles[j].clear();
229 edm::Handle< std::vector<pat::Electron> > ele_handle;
230 iEvent.getByLabel(electron_sources[j], ele_handle);
231 const std::vector<pat::Electron>& pat_electrons = *(ele_handle.product());
232
233 for (unsigned int i = 0; i < pat_electrons.size(); ++i) {
234 pat::Electron pat_ele = pat_electrons[i];
235 Electron ele;
236
237 ele.charge = pat_ele.charge();
238 ele.pt = pat_ele.pt();
239 ele.eta = pat_ele.eta();
240 ele.phi = pat_ele.phi();
241 ele.energy = pat_ele.energy();
242 ele.vertex_x = pat_ele.vertex().x();
243 ele.vertex_y = pat_ele.vertex().y();
244 ele.vertex_z = pat_ele.vertex().z();
245 ele.supercluster_eta = pat_ele.superCluster()->eta();
246 ele.supercluster_phi = pat_ele.superCluster()->phi();
247 ele.dB = pat_ele.dB();
248 //ele.particleIso = pat_ele.particleIso();
249 ele.neutralHadronIso = pat_ele.neutralHadronIso();
250 ele.chargedHadronIso = pat_ele.chargedHadronIso();
251 ele.trackIso = pat_ele.trackIso();
252 ele.photonIso = pat_ele.photonIso();
253 ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
254 ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
255 ele.gsfTrack_px= pat_ele.gsfTrack()->px();
256 ele.gsfTrack_py= pat_ele.gsfTrack()->py();
257 ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
258 ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
259 ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
260 ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
261 ele.passconversionveto = !ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position());
262 ele.dEtaIn=pat_ele.deltaEtaSuperClusterTrackAtVtx();
263 ele.dPhiIn=pat_ele.deltaPhiSuperClusterTrackAtVtx();
264 ele.sigmaIEtaIEta=pat_ele.sigmaIetaIeta();
265 ele.HoverE=pat_ele.hadronicOverEm();
266 ele.fbrem=pat_ele.fbrem();
267 ele.EoverPIn=pat_ele.eSuperClusterOverP();
268 ele.EcalEnergy=pat_ele.ecalEnergy();
269 ele.mvaTrigV0=pat_ele.electronID("mvaTrigV0");
270 ele.mvaNonTrigV0=pat_ele.electronID("mvaNonTrigV0");
271
272 eles[j].push_back(ele);
273 }
274 }
275 }
276
277 // ------------- muons -------------
278 if(doMuons){
279 for(size_t j=0; j< muon_sources.size(); ++j){
280 mus[j].clear();
281
282 edm::Handle< std::vector<pat::Muon> > mu_handle;
283 iEvent.getByLabel(muon_sources[j], mu_handle);
284 const std::vector<pat::Muon>& pat_muons = *(mu_handle.product());
285
286 for (unsigned int i = 0; i <pat_muons.size() ; ++i) {
287 pat::Muon pat_mu = pat_muons[i];
288
289 Muon mu;
290 mu.charge = pat_mu.charge();
291 mu.pt = pat_mu.pt();
292 mu.eta = pat_mu.eta();
293 mu.phi = pat_mu.phi();
294 mu.energy = pat_mu.energy();
295 mu.vertex_x = pat_mu.vertex().x();
296 mu.vertex_y = pat_mu.vertex().y();
297 mu.vertex_z = pat_mu.vertex().z();
298 mu.dB = pat_mu.dB();
299 //mu.particleIso = pat_mu.particleIso();
300 mu.neutralHadronIso = pat_mu.neutralHadronIso();
301 mu.chargedHadronIso = pat_mu.chargedHadronIso();
302 mu.trackIso = pat_mu.trackIso();
303 mu.photonIso = pat_mu.photonIso();
304 mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
305 mu.isGlobalMuon = pat_mu.isGlobalMuon();
306 mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
307 mu.isTrackerMuon = pat_mu.isTrackerMuon();
308 mu.numberOfMatchedStations = pat_mu.numberOfMatchedStations();
309 reco::TrackRef globalTrack = pat_mu.globalTrack();
310 if(!globalTrack.isNull()){
311 mu.globalTrack_chi2 = globalTrack->chi2();
312 mu.globalTrack_ndof = globalTrack->ndof();
313 mu.globalTrack_d0 = globalTrack->d0();
314 mu.globalTrack_d0Error = globalTrack->d0Error();
315 mu.globalTrack_numberOfValidHits = globalTrack->numberOfValidHits();
316 mu.globalTrack_numberOfLostHits = globalTrack->numberOfLostHits();
317 }
318 else{
319 mu.globalTrack_chi2 = 0;
320 mu.globalTrack_ndof = 0;
321 mu.globalTrack_d0 = 0;
322 mu.globalTrack_d0Error = 0;
323 mu.globalTrack_numberOfValidHits = 0;
324 mu.globalTrack_numberOfLostHits = 0;
325 }
326 reco::TrackRef innerTrack = pat_mu.innerTrack();
327 if(!innerTrack.isNull()){
328 mu.innerTrack_chi2 = innerTrack->chi2();
329 mu.innerTrack_ndof = innerTrack->ndof();
330 mu.innerTrack_d0 = innerTrack->d0();
331 mu.innerTrack_d0Error = innerTrack->d0Error();
332 mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
333 mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
334 mu.innerTrack_trackerLayersWithMeasurement = innerTrack->hitPattern().trackerLayersWithMeasurement();
335 mu.innerTrack_numberOfValidPixelHits = innerTrack->hitPattern().numberOfValidPixelHits();
336 }
337 else{
338 mu.innerTrack_chi2 = 0;
339 mu.innerTrack_ndof = 0;
340 mu.innerTrack_d0 = 0;
341 mu.innerTrack_d0Error = 0;
342 mu.innerTrack_numberOfValidHits = 0;
343 mu.innerTrack_numberOfLostHits = 0;
344 mu.innerTrack_trackerLayersWithMeasurement = 0;
345 mu.innerTrack_numberOfValidPixelHits = 0;
346 }
347 reco::TrackRef outerTrack = pat_mu.outerTrack();
348 if(!outerTrack.isNull()){
349 mu.outerTrack_chi2 = outerTrack->chi2();
350 mu.outerTrack_ndof = outerTrack->ndof();
351 mu.outerTrack_d0 = outerTrack->d0();
352 mu.outerTrack_d0Error = outerTrack->d0Error();
353 mu.outerTrack_numberOfValidHits = outerTrack->numberOfValidHits();
354 mu.outerTrack_numberOfLostHits = outerTrack->numberOfLostHits();
355 }
356 else{
357 mu.outerTrack_chi2 = 0;
358 mu.outerTrack_ndof = 0;
359 mu.outerTrack_d0 = 0;
360 mu.outerTrack_d0Error = 0;
361 mu.outerTrack_numberOfValidHits = 0;
362 mu.outerTrack_numberOfLostHits = 0;
363 }
364
365 mus[j].push_back(mu);
366 }
367 }
368 }
369 // ------------- taus -------------
370
371 if(doTaus){
372 for(size_t j=0; j< tau_sources.size(); ++j){
373 taus[j].clear();
374
375
376 edm::Handle< std::vector<pat::Tau> > tau_handle;
377 iEvent.getByLabel(tau_sources[j], tau_handle);
378 const std::vector<pat::Tau>& pat_taus = *(tau_handle.product());
379
380 for (unsigned int i = 0; i < pat_taus.size(); ++i) {
381 pat::Tau pat_tau = pat_taus[i];
382 if(pat_tau.pt() < tau_ptmin) continue;
383 if(fabs(pat_tau.eta()) > tau_etamax) continue;
384
385 Tau tau;
386 tau.charge = pat_tau.charge();
387 tau.pt = pat_tau.pt();
388 tau.eta = pat_tau.eta();
389 tau.phi = pat_tau.phi();
390 tau.energy = pat_tau.energy();
391 tau.decayModeFinding = pat_tau.tauID("decayModeFinding")>0.5;
392 tau.byVLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5;
393 tau.byLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5;
394 tau.byMediumCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5;
395 tau.byTightCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5;
396 tau.againstElectronLoose = pat_tau.tauID("againstElectronLoose")>0.5;
397 tau.againstElectronMedium = pat_tau.tauID("againstElectronMedium")>0.5;
398 tau.againstElectronTight = pat_tau.tauID("againstElectronTight")>0.5;
399 tau.againstElectronMVA = pat_tau.tauID("againstElectronMVA")>0.5;
400 tau.againstMuonLoose = pat_tau.tauID("againstMuonLoose")>0.5;
401 tau.againstMuonMedium = pat_tau.tauID("againstMuonMedium")>0.5;
402 tau.againstMuonTight = pat_tau.tauID("againstMuonTight")>0.5;
403
404 reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
405 if(!leadPFCand.isNull()){
406 tau.leadPFCand_px = leadPFCand->px();
407 tau.leadPFCand_py = leadPFCand->py();
408 tau.leadPFCand_pz = leadPFCand->pz();
409 }
410 else{
411 tau.leadPFCand_px = 0;
412 tau.leadPFCand_py = 0;
413 tau.leadPFCand_pz = 0;
414 }
415 taus[j].push_back(tau);
416 }
417 }
418 }
419
420 // ------------- jets -------------
421 if(doJets){
422 for(size_t j=0; j< jet_sources.size(); ++j){
423
424 jets[j].clear();
425
426 edm::Handle< std::vector<pat::Jet> > jet_handle;
427 iEvent.getByLabel(jet_sources[j], jet_handle);
428 const std::vector<pat::Jet>& pat_jets = *(jet_handle.product());
429
430 for (unsigned int i = 0; i < pat_jets.size(); ++i) {
431 pat::Jet pat_jet = pat_jets[i];
432 if(pat_jet.pt() < jet_ptmin) continue;
433 if(fabs(pat_jet.eta()) > jet_etamax) continue;
434 // std::cout << "available btag discriminators: " << std::endl;
435 // const std::vector<std::pair<std::string, float> > & dis = pat_jets[i].getPairDiscri();
436 // for(size_t k=0; k<dis.size(); ++k){
437 // std::cout << dis[k].first << std::endl;
438 // }
439
440 Jet jet;
441 jet.charge = pat_jet.charge();
442 jet.pt = pat_jet.pt();
443 jet.eta = pat_jet.eta();
444 jet.phi = pat_jet.phi();
445 jet.energy = pat_jet.energy();
446 jet.numberOfDaughters =pat_jet.numberOfDaughters();
447 const reco::TrackRefVector& jettracks = pat_jet.associatedTracks();
448 jet.nTracks = jettracks.size();
449 jet.jetArea = pat_jet.jetArea();
450 jet.pileup = pat_jet.pileup();
451 if(pat_jet.isPFJet()){
452 jet.neutralEmEnergyFraction =pat_jet.neutralEmEnergyFraction();
453 jet.neutralHadronEnergyFraction =pat_jet.neutralHadronEnergyFraction();
454 jet.chargedEmEnergyFraction =pat_jet.chargedEmEnergyFraction();
455 jet.chargedHadronEnergyFraction =pat_jet.chargedHadronEnergyFraction();
456 jet.muonEnergyFraction =pat_jet.muonEnergyFraction();
457 jet.photonEnergyFraction =pat_jet.photonEnergyFraction();
458 jet.chargedMultiplicity =pat_jet.chargedMultiplicity();
459 jet.neutralMultiplicity =pat_jet.neutralMultiplicity();
460 jet.muonMultiplicity =pat_jet.muonMultiplicity();
461 jet.electronMultiplicity =pat_jet.electronMultiplicity();
462 jet.photonMultiplicity =pat_jet.photonMultiplicity();
463 }
464
465 jecUnc->setJetEta(pat_jet.eta());
466 jecUnc->setJetPt(pat_jet.pt());
467 jet.JEC_uncertainty = jecUnc->getUncertainty(true);
468 jet.JEC_factor_raw = pat_jet.jecFactor("Uncorrected");
469
470 jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
471 jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
472 jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
473 jet.btag_combinedSecondaryVertexMVA=pat_jet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
474 jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
475 jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
476
477 const reco::GenJet *genj = pat_jet.genJet();
478 if(genj){
479 jet.genjet_pt = genj->pt();
480 jet.genjet_eta = genj->eta();
481 jet.genjet_phi = genj->phi();
482 jet.genjet_energy = genj->energy();
483 }
484
485 jets[j].push_back(jet);
486 }
487 }
488 }
489
490 // ------------- top jets -------------
491 if(doTopJets){
492 for(size_t j=0; j< topjet_sources.size(); ++j){
493
494 topjets[j].clear();
495
496 edm::Handle<pat::JetCollection> pat_topjets;
497 //edm::Handle<std::vector<reco::Jet> > pat_topjets;
498 iEvent.getByLabel(topjet_sources[j], pat_topjets);
499
500 for (unsigned int i = 0; i < pat_topjets->size(); i++) {
501 const pat::Jet pat_topjet = * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
502 if(pat_topjet.pt() < topjet_ptmin) continue;
503 if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
504
505 TopJet topjet;
506 topjet.charge = pat_topjet.charge();
507 topjet.pt = pat_topjet.pt();
508 topjet.eta = pat_topjet.eta();
509 topjet.phi = pat_topjet.phi();
510 topjet.energy = pat_topjet.energy();
511 topjet.numberOfDaughters =pat_topjet.numberOfDaughters();
512 const reco::TrackRefVector& topjettracks = pat_topjet.associatedTracks();
513 topjet.nTracks = topjettracks.size();
514 topjet.jetArea = pat_topjet.jetArea();
515 topjet.pileup = pat_topjet.pileup();
516 // topjet.neutralEmEnergyFraction =pat_topjet.neutralEmEnergyFraction();
517 // topjet.neutralHadronEnergyFraction =pat_topjet.neutralHadronEnergyFraction();
518 // topjet.chargedEmEnergyFraction =pat_topjet.chargedEmEnergyFraction();
519 // topjet.chargedHadronEnergyFraction =pat_topjet.chargedHadronEnergyFraction();
520 // topjet.muonEnergyFraction =pat_topjet.muonEnergyFraction();
521 // topjet.photonEnergyFraction =pat_topjet.photonEnergyFraction();
522 // topjet.chargedMultiplicity =pat_topjet.chargedMultiplicity();
523 // topjet.neutralMultiplicity =pat_topjet.neutralMultiplicity();
524 // topjet.muonMultiplicity =pat_topjet.muonMultiplicity();
525 // topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
526 // topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
527
528 jecUnc->setJetEta(pat_topjet.eta());
529 jecUnc->setJetPt(pat_topjet.pt());
530 topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
531 topjet.JEC_factor_raw = pat_topjet.jecFactor("Uncorrected");
532
533 topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
534 topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
535 topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
536 topjet.btag_combinedSecondaryVertexMVA=pat_topjet.bDiscriminator("combinedSecondaryVertexMVABJetTags");
537 topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
538 topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
539
540 const reco::GenJet *genj = pat_topjet.genJet();
541 if(genj){
542 topjet.genjet_pt = genj->pt();
543 topjet.genjet_eta = genj->eta();
544 topjet.genjet_phi = genj->phi();
545 topjet.genjet_energy = genj->energy();
546 }
547
548 for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
549 Particle subjet_v4;
550 subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
551 subjet_v4.eta = pat_topjet.daughter(k)->p4().eta();
552 subjet_v4.phi = pat_topjet.daughter(k)->p4().phi();
553 subjet_v4.energy = pat_topjet.daughter(k)->p4().E();
554 topjet.subjets.push_back(subjet_v4);
555 }
556 topjets[j].push_back(topjet);
557 }
558 }
559 }
560
561 // ------------- photons -------------
562 if(doPhotons){
563 for(size_t j=0; j< photon_sources.size(); ++j){
564 phs[j].clear();
565
566 edm::Handle< std::vector<pat::Photon> > photon_handle;
567 iEvent.getByLabel(photon_sources[j], photon_handle);
568 const std::vector<pat::Photon>& pat_photons = *(photon_handle.product());
569
570 for (unsigned int i = 0; i < pat_photons.size(); ++i) {
571 pat::Photon pat_photon = pat_photons[i];
572 Photon ph;
573 ph.charge = 0;
574 ph.pt = pat_photon.pt();
575 ph.eta = pat_photon.eta();
576 ph.phi = pat_photon.phi();
577 ph.energy = pat_photon.energy();
578 ph.vertex_x = pat_photon.vertex().x();
579 ph.vertex_y = pat_photon.vertex().y();
580 ph.vertex_z = pat_photon.vertex().z();
581 ph.supercluster_eta = pat_photon.superCluster()->eta();
582 ph.supercluster_phi = pat_photon.superCluster()->phi();
583 // ph.neutralHadronIso = pat_photon.neutralHadronIso();
584 // ph.chargedHadronIso = pat_photon.chargedHadronIso();
585 ph.trackIso = pat_photon.trackIso();
586 phs[j].push_back(ph);
587 }
588 }
589 }
590
591 // ------------- MET -------------
592 if(doMET){
593 for(size_t j=0; j< met_sources.size(); ++j){
594
595 edm::Handle< std::vector<pat::MET> > met_handle;
596 iEvent.getByLabel(met_sources[j], met_handle);
597 const std::vector<pat::MET>& pat_mets = *(met_handle.product());
598
599 if(pat_mets.size()!=1){
600 std::cout<< "WARNING: number of METs = " << pat_mets.size() <<", should be 1" << std::endl;
601 }
602 else{
603 pat::MET pat_met = pat_mets[0];
604
605 met[j].pt=pat_met.pt();
606 met[j].phi=pat_met.phi();
607 met[j].mEtSig= pat_met.mEtSig();
608 }
609
610 }
611 }
612
613 // ------------- trigger -------------
614 if(doTrigger){
615 edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
616 edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
617 iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
618
619 const edm::Provenance *meta = dummy_TriggerEvent.provenance();
620 std::string nameProcess = meta->processName();
621 edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
622 triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
623
624 edm::Handle<edm::TriggerResults> trigger;
625 iEvent.getByLabel(triggerResultTag, trigger);
626 const edm::TriggerResults& trig = *(trigger.product());
627
628 triggerResults.clear();
629 triggerNames.clear();
630 L1_prescale.clear();
631 HLT_prescale.clear();
632
633 edm::Service<edm::service::TriggerNamesService> tns;
634 std::vector<std::string> triggerNames_all;
635 tns->getTrigPaths(trig,triggerNames_all);
636
637 if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
638 for(unsigned int i=0; i<trig.size(); ++i){
639 std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
640 for(; it!=trigger_prefixes.end(); ++it){
641 if(triggerNames_all[i].substr(0, it->size()) == *it)break;
642 }
643 if(it==trigger_prefixes.end()) continue;
644
645 //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
646 triggerResults.push_back(trig.accept(i));
647 if(newrun) triggerNames.push_back(triggerNames_all[i]);
648 if(isRealData){
649 std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
650 L1_prescale.push_back(pre.first);
651 HLT_prescale.push_back(pre.second);
652 //std::cout << triggerNames_all[i] << " " << pre.first << " " <<pre.second << " " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
653 }
654 }
655 // for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
656 // std::cout << (*iter).first << " " << (*iter).second << std::endl;
657 // }
658 newrun=false;
659 }
660
661 // ------------- generator info -------------
662
663 if(doGenInfo){
664 genInfo.weights.clear();
665 genInfo.binningValues.clear();
666 genps.clear();
667
668 edm::Handle<GenEventInfoProduct> genEventInfoProduct;
669 iEvent.getByLabel("generator", genEventInfoProduct);
670 const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
671
672 genInfo.binningValues = genEventInfo.binningValues();
673 genInfo.weights = genEventInfo.weights();
674 genInfo.alphaQCD = genEventInfo.alphaQCD();
675 genInfo.alphaQED = genEventInfo.alphaQED();
676 genInfo.qScale = genEventInfo.qScale();
677
678 const gen::PdfInfo* pdf = genEventInfo.pdf();
679 if(pdf){
680 genInfo.pdf_id1=pdf->id.first;
681 genInfo.pdf_id2=pdf->id.second;
682 genInfo.pdf_x1=pdf->x.first;
683 genInfo.pdf_x2=pdf->x.second;
684 genInfo.pdf_scalePDF=pdf->scalePDF;
685 genInfo.pdf_xPDF1=pdf->xPDF.first;
686 genInfo.pdf_xPDF2=pdf->xPDF.second;
687 }
688 else{
689 genInfo.pdf_id1=-999;
690 genInfo.pdf_id2=-999;
691 genInfo.pdf_x1=-999;
692 genInfo.pdf_x2=-999;
693 genInfo.pdf_scalePDF=-999;
694 genInfo.pdf_xPDF1=-999;
695 genInfo.pdf_xPDF2=-999;
696 }
697
698 edm::Handle<std::vector<PileupSummaryInfo> > pus;
699 iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
700 genInfo.pileup_NumInteractions_intime=0;
701 genInfo.pileup_NumInteractions_ootbefore=0;
702 genInfo.pileup_NumInteractions_ootafter=0;
703 if(pus.isValid()){
704 genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
705 for(size_t i=0; i<pus->size(); ++i){
706 if(pus->at(i).getBunchCrossing() == 0) // intime pileup
707 genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
708 else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
709 genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
710 }
711 else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
712 genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
713 }
714 }
715 }
716
717 edm::Handle<reco::GenParticleCollection> genPartColl;
718 iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
719 int index=-1;
720 for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
721 index++;
722
723 //write out only top quarks and status 3 particles (works fine only for MadGraph)
724 if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
725 GenParticle genp;
726 genp.charge = iter->charge();;
727 genp.pt = iter->p4().pt();
728 genp.eta = iter->p4().eta();
729 genp.phi = iter->p4().phi();
730 genp.energy = iter->p4().E();
731 genp.index =index;
732 genp.status = iter->status();
733 genp.pdgId = iter->pdgId();
734
735 genp.mother1=-1;
736 genp.mother2=-1;
737 genp.daughter1=-1;
738 genp.daughter2=-1;
739
740 int nm=iter->numberOfMothers();
741 int nd=iter->numberOfDaughters();
742
743 if (nm>0) genp.mother1 = iter->motherRef(0).key();
744 if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
745 if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
746 if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
747
748 genps.push_back(genp);
749 }
750 }
751
752 }
753
754 tr->Fill();
755 if(doLumiInfo)
756 previouslumiblockwasfilled=true;
757 }
758
759
760 // ------------ method called once each job just before starting event loop ------------
761 void
762 NtupleWriter::beginJob()
763 {
764 if(doLumiInfo){
765 totalRecLumi=0;
766 totalDelLumi=0;
767 previouslumiblockwasfilled=false;
768 }
769 }
770
771 // ------------ method called once each job just after ending the event loop ------------
772 void
773 NtupleWriter::endJob()
774 {
775 outfile->cd();
776 tr->Write();
777 outfile->Close();
778 }
779
780 // ------------ method called when starting to processes a run ------------
781 void
782 NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup)
783 {
784 if(doTrigger){
785 bool setup_changed = false;
786 hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
787 newrun=true;
788 }
789
790 if(doJets || doTopJets){
791 edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
792 iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
793 JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
794 jecUnc = new JetCorrectionUncertainty(JetCorPar);
795 }
796 }
797
798 // ------------ method called when ending the processing of a run ------------
799 void
800 NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
801 {
802 if(doLumiInfo)
803 std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
804 }
805
806 // ------------ method called when starting to processes a luminosity block ------------
807 void
808 NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
809 {
810 if(doLumiInfo){
811 edm::Handle<LumiSummary> l;
812 lumi.getByLabel("lumiProducer", l);
813
814 //add lumi of lumi blocks without any event to next lumiblock
815 if(previouslumiblockwasfilled){
816 intgRecLumi=0;
817 intgDelLumi=0;
818 }
819 previouslumiblockwasfilled=false;
820
821 if (l.isValid()){;
822 intgRecLumi+=l->intgRecLumi()*6.37;
823 intgDelLumi+=l->intgDelLumi()*6.37;
824 totalRecLumi+=l->intgRecLumi()*6.37;
825 totalDelLumi+=l->intgDelLumi()*6.37;
826 }
827 //std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<" " << l->intgDelLumi()*6.37<<std::endl;
828 //std::cout << "summed: "<< intgRecLumi << " " << intgDelLumi << std::endl;
829 }
830 }
831
832 // ------------ method called when ending the processing of a luminosity block ------------
833 void
834 NtupleWriter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
835 {
836 }
837
838 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
839 void
840 NtupleWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
841 //The following says we do not know what parameters are allowed so do no validation
842 // Please change this to state exactly what you do use, even if it is no parameters
843 edm::ParameterSetDescription desc;
844 desc.setUnknown();
845 descriptions.addDefault(desc);
846 }
847
848 //define this as a plug-in
849 DEFINE_FWK_MODULE(NtupleWriter);