ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
Revision: 1.9
Committed: Mon Apr 16 15:42:09 2012 UTC (13 years ago) by peiffer
Content type: text/plain
Branch: MAIN
Changes since 1.8: +14 -2 lines
Log Message:
tau id and ttbsm cfg

File Contents

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